package stats;
import java.util.*;

public class ProbLib {
	
	static Random generator = new Random();
	
	public static double meanMax(double mean1,double stdev1, double mean2,double stdev2,double ro)
	{
		double var1 = Math.pow(stdev1,2);
		double var2 = Math.pow(stdev2,2);
		double teta = Math.sqrt(Math.pow(var1, 2) + Math.pow(var2, 2) - 2*ro*var1*var2);
		//double teta = Math.sqrt(stdev1 + stdev2 - 2*ro*Math.sqrt(stdev1*stdev2));
		return mean1 * Phi((mean1-mean2)/teta) + mean2 * Phi((mean2-mean1)/teta) + teta * phi((mean1-mean2)/teta);
	}

	public static double meanMaxCheck(double mean1,double stdev1, double mean2,double stdev2,double ro)
	{
		int trials = 100000;
		double moment=0;
		for (int i=0;i<trials;i++)
		{
			double rnd1 = generateGaussian(mean1,stdev1);
			double rnd2 = generateGaussian(mean2,stdev2);
			if (rnd1>rnd2)
			{
				moment+=rnd1;
			}
			else
			{
				moment+=rnd2;
			}
		}
		return moment/trials;
	}
	
	public static double stdevMax(double mean1,double stdev1, double mean2,double stdev2,double ro)
	{
		double var1 = Math.pow(stdev1,2);
		double var2 = Math.pow(stdev2,2);
//		double teta = Math.sqrt(Math.pow(var1, 2) + Math.pow(var2, 2) - 2*ro*var1*var2);
//		return Math.sqrt((Math.pow(var1,2) + Math.pow(mean1,2)) * Phi((mean1-mean2)/teta) + (Math.pow(var2,2) + Math.pow(mean2,2)) * Phi((mean2-mean1)/teta) + (mean1+mean2) * teta * phi((mean1-mean2)/teta));
		double teta = Math.sqrt(Math.pow(stdev1, 2) + Math.pow(stdev2, 2) - 2*ro*stdev1*stdev2);
		return (Math.pow(stdev1,2) + Math.pow(mean1,2)) * Phi((mean1-mean2)/teta) + (Math.pow(stdev2,2) + Math.pow(mean2,2)) * Phi((mean2-mean1)/teta) + (mean1+mean2) * teta * phi((mean1-mean2)/teta);
	}	
	
	public static double stdevMaxCheck(double mean1,double stdev1, double mean2,double stdev2,double ro)
	{
		int trials = 100000;
		
		double rndMax[] = new double[trials];
		
		double mean=0;
		for (int i=0;i<trials;i++)
		{
			double rnd1 = generateGaussian(mean1,stdev1);
			double rnd2 = generateGaussian(mean2,stdev2);
			if (rnd1>rnd2)
			{
				rndMax[i]=rnd1;
			}
			else
			{
				rndMax[i]=rnd2;
			}
			mean+=rndMax[i];
		}
		mean = mean/trials;

		double stdev=0;
		for (int i=0;i<trials;i++)
		{
			stdev += Math.pow(rndMax[i]-mean, 2);
		}
		
		return Math.sqrt(stdev/(trials-1));
	}

	public static double partialPositiveFirstMoment(double mean,double stdev)
	{
		return Math.pow(stdev, 2)*phi(mean/stdev) + mean * (1 - Phi(-mean/stdev));
	}	

	public static double partialPositiveFirstMomentCheck(double mean,double stdev)
	{
		int trials = 100000;
		double moment=0;
		for (int i=0;i<trials;i++)
		{
			double rnd = generateGaussian(mean,stdev);
			if (rnd>0)
			{
				moment+=rnd;
			}
		}
		return moment/trials;
	}

	public static double partialNegativeFirstMoment(double mean,double stdev)
	{
		return - Math.pow(stdev, 2)*phi(mean/stdev) + mean * Phi(-mean/stdev);
	}		

	public static double partialNegativeFirstMomentCheck(double mean,double stdev)
	{
		int trials = 100000;
		double moment=0;
		for (int i=0;i<trials;i++)
		{
			double rnd = generateGaussian(mean,stdev);
			if (rnd<0)
			{
				moment+=rnd;
			}
		}
		return moment/trials;
	}

	public static double phi(double x)
	{
		return 1 / Math.sqrt(2*Math.PI) * Math.pow(Math.E, - Math.pow(x,2)/2);
	}

	public static double Phi(double x)
	{
		return CDF_Normal.normp(x); 
	}

	public static double PhiCheck(double x)
	{
		int trials = 10000;
		int successes=0;
		for (int i=0;i<trials;i++)
		{
			if (generateGaussian(0,1)<x)
			{
				successes++;
			}
		}
		return (double)successes/trials;
	}
	
	public static double generateGaussian(double mean, double stdev)
	{
		return mean + Math.pow(stdev,2) * generator.nextGaussian();
	}
}
