Copyright (C) 2005-2009 Chad Brassil; Last modified 10 June 2009

Set-up

Execute this section first

(* These commands turn off some spelling warning messages *)
Off[General::"spell1"];
Off[General::"spell"];
Off[Graphics::gptn];

(****** OTHER FUNCTIONS ******)

Clear[CPlot];
CPlot[x_List,y_List,opt__]:=ListPlot[Transpose[{x,y}],opt,PlotStyle->PointSize[0.02`],Frame->True,FrameLabel->{If[Length[Position[Transpose[data],x]]≠0,labels[[Position[Transpose[data],x][[1,1]]]],""],If[Length[Position[Transpose[data],y]]≠0,labels[[Position[Transpose[data],y][[1,1]]]],""]},BaseStyle->{FontFamily->"Swiss721BT-Roman",FontSize->12},PlotRegion->{{0,1},{0,0.99`}},ImageSize->350]
CPlot[x_List,y_List]:=ListPlot[Transpose[{x,y}],PlotStyle->PointSize[0.02`],Frame->True,FrameLabel->{If[Length[Position[Transpose[data],x]]≠0,labels[[Position[Transpose[data],x][[1,1]]]],""],If[Length[Position[Transpose[data],y]]≠0,labels[[Position[Transpose[data],y][[1,1]]]],""]},BaseStyle->{FontFamily->"Swiss721BT-Roman",FontSize->12},PlotRegion->{{0,1},{0,0.99`}},ImageSize->350]
CPlot::"usage"="CPlot[ x, y, opt] produces a scatter plot given two column labels from the data set.  Standard ListPlot options can also be included";


CPlot3D[x_List,y_List,z_List,opt__]:=Show[Graphics3D[{PointSize[0.02`],Point/@Transpose[{x,y,z}]}],opt,BoxRatios->{1.5`,1.5`,1},Axes->True,AxesLabel->{If[Length[Position[Transpose[data],x]]≠0,labels[[Position[Transpose[data],x][[1,1]]]],""],If[Length[Position[Transpose[data],y]]≠0,labels[[Position[Transpose[data],y][[1,1]]]],""],If[Length[Position[Transpose[data],z]]≠0,labels[[Position[Transpose[data],z][[1,1]]]],""]},AxesEdge->{{-1,-1},{1,-1},{-1,-1}},BaseStyle->{FontFamily->"Swiss721BT-Roman",FontSize->12},ImageSize->350]
CPlot3D[x_List,y_List,z_List]:=Show[Graphics3D[{PointSize[0.02`],Point/@Transpose[{x,y,z}]}],BoxRatios->{1.5`,1.5`,1},Axes->True,AxesLabel->{If[Length[Position[Transpose[data],x]]≠0,labels[[Position[Transpose[data],x][[1,1]]]],""],If[Length[Position[Transpose[data],y]]≠0,labels[[Position[Transpose[data],y][[1,1]]]],""],If[Length[Position[Transpose[data],z]]≠0,labels[[Position[Transpose[data],z][[1,1]]]],""]},AxesEdge->{{-1,-1},{1,-1},{-1,-1}},BaseStyle->{FontFamily->"Swiss721BT-Roman",FontSize->12},ImageSize->350]
CPlot3D::"usage"="CPlot3D[ x, y, z, opt] produces a scatter plot given three column labels from the data set.  Standard Graphics3D options can also be included";


Categorical[parameter_,categoryLabel_List,opt___]:=Block[{i,tempCategoryLabel},If[Length[Dimensions[categoryLabel]]>1,tempCategoryLabel=StringJoin/@Transpose[categoryLabel],tempCategoryLabel=categoryLabel;];Flatten[Table[(Subscript[parameter, Symbol[tempCategoryLabel[[i]]]]),{i,Length[tempCategoryLabel]}]/.{opt}]];
Categorical::"usage"="Categorical[ parameter, categoryLabel, optionalReplacementRules ] generates category specific versions of parameter based on the categories listed in categoryLabel.  To examine interactions enter a list for categoryLabel, for example {factorA,factorB}.  Optional replacement rules (optionalReplacementRules) can be used to combine categories; for example {a->ab,b->ab}";


CategoricalMean[dataColumn_List,categoryLabel_List,opt___]:=Block[{i,tempCategoryLabel,categoricalList,meansPerCat},tempCategoryLabel=categoryLabel/.{opt};If[Length[Dimensions[tempCategoryLabel]]>1,tempCategoryLabel=StringJoin/@Transpose[tempCategoryLabel]];categoricalList=Symbol/@Column[Reverse[Sort[Tally[tempCategoryLabel]],2],2];meansPerCat=(Mean[SelectCases[dataColumn,tempCategoryLabel,#1]]&)/@categoricalList;Thread[categoricalList->meansPerCat]];
CategoricalMean::"usage"="CategoricalMean[ dataColumn, categoryLabel, optionalReplacementRules ] generates category specific Means of a dataColumn based on the categories listed in categoryLabel.  To examine interactions enter a list for categoryLabel, for example {factorA,factorB}.  Optional replacement rules (optionalReplacementRules) can be used to combine categories; for example {a->ab,b->ab}";


CategoricalSymbols[categoryLabel_List]:=Symbol/@Column[Reverse[Sort[Tally[categoryLabel]],2],2];
CategoricalSymbols::"usage"="CategoricalSymbols[ categoryLabel ] generates a list of categories in categoryLabel.  Output can be used as base for a substitution list along with CategoricalMean output.";


SelectCases[dataColumn_List,categoryColumn_List,equalToCriteria_]:=If[Length[Dimensions[categoryColumn]]==1,Extract[dataColumn,Position[(TrueQ[Or@@Thread[#1==Flatten[{equalToCriteria}]]]&)/@ToExpression[categoryColumn],True]],
Extract[dataColumn,Position[(TrueQ[Or@@Thread[#1==Flatten[{Symbol[StringJoin[{equalToCriteria}]]}]]]&)/@ToExpression[StringJoin/@(categoryColumn\[Transpose])],True]]
];
(* This could still use work to make it robust *)
SelectCases::"usage"="SelectCases[ dataColumn, categoryColumn, equalToCriteria] produces a list of data values from dataColumn where categoryColumn == equalToCriteria.  To select multiple criteria, use a list of equalToCriteria, for example {a,b}.  To select multiple categories, use a list in for categoryColumn and equaltoCriteria";


ProcessData[dataInput_]:=Module[{},If[NumericQ[dataInput[[1,1]]]==False,labels=dataInput[[1]];data=Drop[dataInput,1];];Clear@@labels;Do[Evaluate[Symbol[labels[[i]]]]=data[[All,i]],{i,Length[labels]}]];
ProcessData::"usage"="ProcessData[ data] takes a data matrix with text column labels (in quotes) and extracts column headings.  Column headings are output into labels with the remaining data being left in data";



ProcessVars[ll_]:=Tally[Select[Flatten[ll/.{sfPoisson->({##1}&),sfNegBinom->({##1}&),sfGenPoisson->({##1}&),sfBinomial->({##1}&),sfNormal->({##1}&),sfBimodalNormal->({##1}&),sfLogNormal->({##1}&),sfBetaMean->({##1}&),sfBetaMeanShape->({##}&),sfBetaBinomial->({##1}&),Plus->({##1}&),Times->({##1}&),Power->({##1}&),Log->({##1}&),"NA"->0}],!NumericQ[#1]&]][[All,1]];




ProcessConstraints[ll_]:=Block[{output},Switch[Head[First[DownValues[ll][[1,2,2]]]],
sfNormal,If[TrueQ[Head[DownValues[ll][[1,2,2,3]]]==List],And@@Thread[Column[Tally[DownValues[ll][[1,2,2,3]]],2]>0],DownValues[ll][[1,2,2,3]]>0],
sfBinomial,If[TrueQ[Head[DownValues[ll][[1,2,2,3]]]==List],And@@Thread[0≤Column[Tally[DownValues[ll][[1,2,2,3]]],2]≤1],0≤DownValues[ll][[1,2,2,3]]≤1],
sfBetaBinomial,If[TrueQ[Head[DownValues[ll][[1,2,2,3]]]==List],output=And@@Thread[0≤Column[Tally[DownValues[ll][[1,2,2,3]]],2]≤1],output=0≤DownValues[ll][[1,2,2,3]]≤1];If[TrueQ[Head[DownValues[ll][[1,2,2,4]]]==List],output=output&&And@@Thread[0<Column[Tally[DownValues[ll][[1,2,2,4]]],2]],output=output&&0<DownValues[ll][[1,2,2,4]]];output,
_,True]];



sfPunishmentConstant=10000.`;




sfPoisson[z_,λ_?NumberQ]:=If[z=="NA",0,If[λ>0,λ-z Log[λ]+Log[z!],sfPunishmentConstant]];
SetAttributes[sfPoisson,Listable];
sfPoisson::"usage"="sfPoisson[ data, λ] is the support function for the Possion distribution where data is an integer value ≥ 0 and λ is the mean";


Clear[sf];
sfNegBinom[z_,m_?NumberQ,k_?NumberQ]:=If[z=="NA",0,If[m>0&&k>0,-Log[Gamma[k+z]]+Log[z!]+Log[Gamma[k]]-k Log[k]+k Log[k+m]-z Log[m]+z Log[k+m],sfPunishmentConstant]];
SetAttributes[sfNegBinom,Listable];
sfNegBinom::"usage"="sfNegBinom[ data, m, k] is the support function for the Negative Binomial distribution where data is an integer value ≥ 0 and m is the mean and k is the shape parameter";


sfGenPoisson[z_?NumberQ,λ1_?NumberQ,λ2_?NumberQ]:=If[z=="NA",0,If[λ1+z λ2>0,-Log[λ1]-(z-1) Log[λ1+z λ2]+(λ1+z λ2)+Log[z!],sfPunishmentConstant]];
SetAttributes[sfGenPoisson,Listable];
sfGenPoisson::"usage"="sfGenPoisson[ data, λ1, λ1] is the support function for the Generalized Possion distribution where data is an integer value ≥ 0 and λ1 is the rate of moving from 0 to 1, and λ2 is the slope of how that rate changes with subsequent categories";



Unprotect[Power];
((x_)^0):=1;

Unprotect[Log];
Log[0.`]:=-∞;
sfBinomial[z_?NumberQ,n_?NumberQ,p_?NumberQ]:=If[z=="NA",0,-Log[n!]+Log[z!]+Log[(n-z)!]-Log[(p^z) (((1.`-p))^(n-z))]/.{∞->100000000000000}];

SetAttributes[sfBinomial,Listable];
sfBinomial::"usage"="sfBinomial[ success, trials, p] is the support function for the Binomial distribution where sucess is an integer value representing the number of successes in integer trails.  p is the probability of success for each trial";





sfBetaBinomial[z_?NumberQ,n_?NumberQ,m_?NumberQ,s_?NumberQ]:=If[z=="NA",0,If[0<m<1&&s>0,Log[(n-z)!]+Log[z!]-Log[n!]+Log[Beta[s m,s (1-m)]]-Log[Beta[s m+z,n+s (1-m)-z]]/.{∞->100000000000000},sfPunishmentConstant]];

SetAttributes[sfBetaBinomial,Listable];
sfBetaBinomial::"usage"="sfBetaBinomial[ success, trials, m, s] is the support function for the BetaBinomial distribution where sucess is an integer value representing the number of successes in integer trails.  Probability of sucess is Beta distributed with a mean probability m and a dispersion parameter s.";


sfLogNormal[z_,m_?NumberQ,σ_?NumberQ]:=If[z=="NA",0,If[σ>0&&m>0,Log[z (Sqrt[2 π (σ^2)])]+((((Log[z]-Log[m]))^2)/(2 (σ^2))),sfPunishmentConstant]];
SetAttributes[sfLogNormal,Listable];
sfLogNormal::"usage"="sfLogNormal[ data, m, σ] is the support function for the Log Normal distribution where data is a positive real number > 0.  m is the median and σ is the standard deviation of the ln(data)";


sfNormal[z_,m_,σ_?NumberQ]:=If[z=="NA",0,Log[(Sqrt[2 π (σ^2)])]+((((z-m))^2)/(2 (σ^2)))];
SetAttributes[sfNormal,Listable];
sfNormal::"usage"="sfNormal[ data, μ, σ] is the support function for the Normal distribution where data is any real number.  μ is the mean and σ is the standard deviation";


sfBimodalNormal[z_,m1_,σ1_?NumberQ,m2_,σ2_?NumberQ]:=If[z=="NA",0,Log[2 (Sqrt[2 π])]-Log[((Exp[-((((z-m1))^2)/(2 ((σ1)^2)))])/(Sqrt[((σ1)^2)]))+((Exp[-((((z-m2))^2)/(2 ((σ2)^2)))])/(Sqrt[((σ2)^2)]))]];
SetAttributes[sfBimodalNormal,Listable];
sfBimodalNormal::"usage"="sfBimodalNormal[ data, μ1, σ1, μ2, σ2] is the support function for the Bimodal Normal distribution where data is any real number.  μ1 is the mean and σ1 is the standard deviation of one mode and μ2 is the mean and σ2 is the standard deviation of the other mode.  Note, this formulation assumes equal weighting of the 2 modes.";


sfBetaMean[z_,μ_?NumberQ,σ_?NumberQ]:=If[z=="NA",0,If[(σ^2)<μ (1-μ)&&0≤μ≤1&&σ>0,(1-(((-1+μ) (-μ+(μ^2)+(σ^2)))/(σ^2))) Log[1-z]+(((-(μ^2)+(μ^3)+(σ^2)+μ (σ^2)) Log[z])/(σ^2))+Log[Beta[-((μ (-μ+(μ^2)+(σ^2)))/(σ^2)),(((-1+μ) (-μ+(μ^2)+(σ^2)))/(σ^2))]],sfPunishmentConstant]];
SetAttributes[sfBetaMean,Listable];
sfBetaMean::"usage"="sfBetaMean[ data, μ, σ] is the support function for the Beta distribution which has been transformed to explicity model the mean (μ) and standard deviation (σ).  Data is a real number between 0 and 1";

(* Beta distribution with mean and shape parameter*)
sfBetaMeanShape[z_,μ_?NumberQ,s_?NumberQ]:=If[z=="NA",0,If[0≤μ≤1&&s>0,
Log[Beta[s μ,s (1-μ)]]-(s μ-1)Log[z]-(s (1-μ)-1)Log[1-z] ,sfPunishmentConstant]];
SetAttributes[sfBetaMeanShape,Listable];
sfBetaMeanShape::usage="sfBetaMeanShape[ data, μ, s] is the support function for the Beta distribution which has been transformed to explicity model the mean (μ) and a shape parameter (s).  Data is a real number between 0 and 1";

CStringQ[column_List]:=And@@Table[TrueQ[Head[column[[i]]]==String],{i,Length[column]}];


ProcessCategorical[categoryHeading_]:=Block[{i,categoryNum,categoricalValues},If[!CStringQ[categoryHeading],Clear[dummyList];Return[False]];categoricalValues=Column[Reverse[Sort[Tally[categoryHeading]],2],2];dummyList={};Do[Clear[Evaluate[categoricalValues[[categoryNum]]<>"$"]];Evaluate[ToExpression[categoricalValues[[categoryNum]]<>"$"]]=Table[If[categoryHeading[[i]]==categoricalValues[[categoryNum]],1,0],{i,Length[categoryHeading]}];AppendTo[dummyList,categoricalValues[[categoryNum]]<>"$"];,{categoryNum,Length[categoricalValues]}];True];

Reference: some basic distributions

An excellent reference for probability distributions is: Evans, Merran, Nicholas Hastings, and Brian Peacock.  2000.  Statistical Distributions, 3rd Edition.  John Wiley & Sons. New York.

Continuous

Expected, or mean
	E[Z]=∫z f(z)dz	
Variance
	VAR(Z)=E[(((Z-E[Z]))^2)]=∫(((z-E[z]))^2)f(z)dz 	

For continuous variables, the probability density function is the probability of the value z given the parameters

Uniform

Uniform Distribution:
a uniform distribution on the interval [a,b] where a < b
probability density function:
	f(z)=(1/(b-a)), where a ≤ z ≤ b
mean: ((b+a)/2)
variance: ((((b-a))^2)/12)

r

dunif(z,min = a, max = b)

# generate 1 uniform random number over the range 1.4 to 2.4
runif(1,1.4,2.4)

mathematica

f[z_,a_,b_]:=(1/(b-a))

RandomReal[]

(* random number of specified range *)

RandomReal[{3,6}]

Normal

Normal Distribution:
has two parameters: m is the location parameter, in this case the mean; σ is the scale parameter, in this case the standard deviation
probability density function:
	f(z)=(1/(Sqrt[2 π (σ^2)]))exp(-((((z-m))^2)/(2 (σ^2))))
mean: m
variance: (σ^2)
note:  plotting f(z) will give you the familiar bell curve

An aside: If one transforms data and then peforms classical statistics assuming a normal distribution, the AIC values cannot be directly compared.  In essence one has changed the data, which is forbidden.  However, a simple transformation of the Likelihood values can render the AIC values comparable.  See <http://www.unc.edu/courses/2006spring/ecol/145/001/docs/lectures/lecture18.htm> for a full derivation.  I quote:

In general, if Z=T(y) is some transformation of the response variable such that T(y) has probability density function f(z).  The likelihood as a function of y is given by the following expression.

L(y)=Π f(T((Subscript[y, i])))(d/(d (Subscript[y, i])))T((Subscript[y, i]))

Remember this is the likelihood, not the negative log likelihood.  Working this through a log-transformed data set means that you simple add Σ 2 ln((Subscript[y, i])) to the AIC value from the transformed analysis, and you will have an AIC that can be compared to analysis on the untransformed data.

r

probability density function

dnorm(z,mean,sd)
#plot the pdf
curve(dnorm(x,4,1),0,7)

random number generation, where n is the number random numbers to generation.  In other words, n = 1 will create 1 random number.

rnorm(n,mean,sd)

mathematica

f[z_,m_,σ_]:=(1/(Sqrt[2π (σ^2)]))Exp[-((((z-m))^2)/(2(σ^2)))]

Plot[f[z,4,1],{z,0,7}]

Show that pdf integrates to 1

Integrate[f[z,1,1],{z,-∞,∞}]

<<Statistics`ContinuousDistributions`
(* normal random number with mean of 3 and standard deviation 2 *)
Random[NormalDistribution[3,2]]

Lognormal

Lognormal Distribution:
has two parameters: m is the location parameter; σ is the shape parameter
probability density function:
	f(z)=(1/(z(Sqrt[2 π (σ^2)])))exp(-((((ln(z)-log(m)))^2)/(2 (σ^2)))), where z > 0, m >0, and σ >0
mean: m exp(((σ^2)/2))
median: m
(σ^2) is the variance of the ln(z)
variance: (m^2)exp((σ^2))(exp((σ^2))-1)
note:  This describes a variable in which the logarithm of the variable is a normal distribution.  Whereas the normal distribution is additive, the lognormal distribution is multiplicative.

r

dlnorm(z,meanlog = m, sdlog = σ)

R uses a slightly different form of the lognormal model.  Use the following form so that m is defined in the same way as described in the f(z) equation above.

dlnorm(z,log( m), σ)

#generate random number
rlnorm(1,meanlog = 1, sdlog = 1.2)

mathematica

f[z_,m_,σ_]:=(1/(z(Sqrt[2π (σ^2)])))Exp[-((((Log[z]-Log[m]))^2)/(2(σ^2)))]

Plot[f[z,1.7`,1],{z,0,20}]

Show that pdf integrates to 1

Integrate[f[z,1,1],{z,0,∞}]

Derive expected value

Integrate[z f[z,m,σ],{z,0,∞},Assumptions->{σ>0,m>0,σ∈Reals,m∈Reals}]

Derive variance

Integrate[(((z-Integrate[z f[z,m,σ],{z,0,∞},Assumptions->{σ>0,m>0,σ∈Reals,m∈Reals}]))^2)f[z,m,σ],{z,0,∞},Assumptions->{σ>0,m>0,σ∈Reals,m∈Reals}]

<<Statistics`ContinuousDistributions`
(* lognormal random number with the first number being the Log[m] and the second number being σ.  In the case below, m = 3 and σ = 2 *)
Random[LogNormalDistribution[Log[3],2]]

Exponential

Exponential Distribution:
has three parameters: m is the location parameter; σ is the shape parameter
probability density function:
	f(z)=a Exp[-a z], where z > 0
mean: (1/a)
variance: (1/(a^2))
note:  Continuous distribution analog of Poisson.  Can be used for dispersal distance.  Widely used in Baysian Analysis

r

dexp(z,rate = a)

#generate random number
rexp(1,rate = a)

mathematica

f[z_,a_]:=a Exp[-a z]

Plot[f[z,0.5`],{z,0,15},PlotRange->All]

<<Statistics`ContinuousDistributions`
(* exponential random number with a = 0.5 *)
Random[ExponentialDistribution[0.5]]

Gamma

Gamma Distribution:
has two parameters: a is the rate parameter; n is the shape parameter
probability density function:
	f(z)=((a^n)/(Γ(n))) Exp[-a z](z^(n-1)), where z > 0
mean: (n/a)
variance: (n/(a^2))
note: can have a long tail.  When n = 1, the gamma distribution becomes the exponential distribution.  When n = degrees of freedom / 2 and a = 2, the gamma distribution becomes the chi-square distribution.

r

dgamma(z,shape = n ,rate = a)

#generate random number
rgamma(1,shape = n,rate = a)

mathematica

f[z_,a_,n_]:=((a^n)/(Gamma[n])) Exp[-a z](z^(n-1))

Plot[f[z,1,1],{z,0,12},PlotRange->All]

Plot[f[z,1,2],{z,0,12},PlotRange->All]

Plot[f[z,0.5`,2],{z,0,12},PlotRange->All]

Derive expected value

Integrate[z f[z,a,n],{z,0,∞},Assumptions->{a>0,n>0,a∈Reals,n∈Reals}]

Derive variance

Integrate[(((z-Integrate[z f[z,a,n],{z,0,∞},Assumptions->{a>0,n>0,a∈Reals,n∈Reals}]))^2)f[z,a,n],{z,0,∞},Assumptions->{a>0,n>0,a∈Reals,n∈Reals}]

<<Statistics`ContinuousDistributions`
(* gamma random number with n = 2 and a = 0.5 *)
Random[GammaDistribution[2.0 ,1/0.5]]

Chi-Square

Note: this section is a bit incomplete

Chi-Square Distribution:  A squared normal distribution is a chi-square distribution.

probability density function for 1 degree of freedom:
	f(z)=(1/(Sqrt[2 π (σ^2)z]))exp(-(z/(2 (σ^2))))
	
probability density function for n degrees of freedom:
	f(z)=(1/((σ^n)(2^((n/2)))Γ((n/2))))(z^((n/2)-1)) exp(-(z/(2 (σ^2)))), for z≥0	
mean: n (σ^2)
variance: 2 n (σ^4)

notes:
-  n can also be thought of as a shape parameter in the context of a general distribution
-  often assume unit variance, σ = 1 
-  Γ((1/2))=(Sqrt[π])

r

dchisq()

mathematica

f[z_,σ_]:=(1/(Sqrt[2π (σ^2)z]))Exp[-(z/(2(σ^2)))]

Plot[f[z,2.1`],{z,0,7}]

Show that pdf integrates to 1

Integrate[f[z,2.1],{z,0,∞}]

<<Statistics`ContinuousDistributions`
(* Chi-Square distributed random number with 1 degree of freedom *)
Random[ChiSquareDistribution[1]]

Add non-centrality, or location, parameter with a single degree of freedom, although no scale parameter
mean: 1 + μ

Two forms are given below

f[z_,μ_]:=Exp[-((z+μ)/2)](1/2)((((z/μ)))^((-1)/4))BesselI[-1/2,(Sqrt[μ z])]

f[z_,μ_]:=Hypergeometric0F1[(1/2),(1/4)z μ](1/(Sqrt[2π z ]))Exp[-((z+μ )/2)]

NIntegrate[ f[z,2.0],{z,0,100000}]

NIntegrate[z f[z,2.0],{z,0,100000}]

Integrate[z f[z,μ],{z,0,∞},Assumptions->{μ>0,μ∈Reals}]

Hypergeometric0F1[(1/2),(1/4)z μ]

(1/2)(Exp[(Sqrt[z μ])]+Exp[-(Sqrt[z μ])])(1/(Sqrt[2π z ]))Exp[-((z+μ )/2)]

(Exp[(Sqrt[z μ])]+Exp[-(Sqrt[z μ])])(1/(2(Sqrt[2π z ])))Exp[-((z+μ )/2)]

(1/(2(Sqrt[2π z ])))Exp[(Sqrt[z μ])-((z+μ )/2)]+(1/(2(Sqrt[2π z ])))Exp[-(Sqrt[z μ])-((z+μ )/2)]

Bimodal Normal

Normal Distribution:
has two parameters: m is the location parameter, in this case the mean; σ is the scale parameter, in this case the standard deviation
probability density function:
	f(z)=q(1/(Sqrt[2 π ((σ1)^2)]))exp(-((((z-m1))^2)/(2 ((σ1)^2))))+(1-q)(1/(Sqrt[2 π ((σ2)^2)]))exp(-((((z-m2))^2)/(2 ((σ2)^2))))
mean of each mode: m1,m2
variance of each mode: ((σ1)^2),((σ2)^2)
note:  plotting f(z) will give you the familiar bell curve 

Here I extend that to a bimodal normal distribution.

r

#assuming equal weighting to each mode
(dnorm(z,m1,s1)+dnorm(z,m2,s2))/2

# weighting c for modes
(1/(1+exp(-c)))*dnorm(z,m1,s1)+(1/(1+exp(c)))dnorm(z,m2,s2)

mathematica

f[z_,m1_,σ1_,m2_,σ2_]:=(1/(2(Sqrt[2π ((σ1)^2)])))Exp[-((((z-m1))^2)/(2((σ1)^2)))]+(1/(2(Sqrt[2π ((σ2)^2)])))Exp[-((((z-m2))^2)/(2((σ2)^2)))]

Plot[f[z,5,1,2,0.7`],{z,0,8}]

Show that pdf integrates to 1

Integrate[f[z,5,1,2,0.7],{z,-∞,∞}]

Therefore the support function is

(* bimodal normal distribution *)
sf[z_,m1_,σ1_,m2_,σ2_]:=Log[2(Sqrt[2 π ])]-Log[(1/(Sqrt[ ((σ1)^2)]))Exp[-((((z-m1))^2)/(2((σ1)^2)))]+(1/(Sqrt[((σ2)^2)]))Exp[-((((z-m2))^2)/(2((σ2)^2)))]];

bimodalRandom[m1_,σ1_,m2_,σ2_]:=Block[{indef,x},indef[x_Real]:=NIntegrate[bimodalPDF[z,m1,σ1,m2,σ2],{z,-∞,x}];x/.FindRoot[RandomReal[]==indef[x],{x,(1/2) (m1+m2) 0.5`,(1/2) (m1+m2) 1.`}]]

Perhaps I should not assume an equal distribution of each mode, in which case I add the parameter c which is between [0,1] and gives the percent of each mode in generating the overall distribution

f[z_,m1_,σ1_,m2_,σ2_,c_]:=(c/(Sqrt[2π ((σ1)^2)]))Exp[-((((z-m1))^2)/(2((σ1)^2)))]+((1-c)/(Sqrt[2π ((σ2)^2)]))Exp[-((((z-m2))^2)/(2((σ2)^2)))]

Plot[f[z,5,1,2,0.7`,0.3`],{z,0,8}]

Show that pdf integrates to 1

Integrate[f[z,5,1,2,0.7,0.3],{z,-∞,∞}]

(* bimodal normal distribution *)
sf[z_,m1_,σ1_,m2_,σ2_,c_]:=Log[(Sqrt[2 π ])]-Log[(c/(Sqrt[ ((σ1)^2)]))Exp[-((((z-m1))^2)/(2((σ1)^2)))]+((1-c)/(Sqrt[((σ2)^2)]))Exp[-((((z-m2))^2)/(2((σ2)^2)))]];

Beta

Beta Distribution:
This distribution is good for modeling fractions because it is constrained to the domain between 0 and 1.  It has two parameters: both p and q are shape parameters with p>0 and q>0, but z is restrict to the range (0,1)

probability density function:
	f(z)=((((z^(p-1))(1-z))^(q-1))/(B(p,q))), B is the beta function, B(p,q) = (∫_0^1)(((u^(p-1))(1-u))^(q-1))du)
	f(z)=((((z^(p-1))(1-z))^(q-1))/((∫_0^1)(((u^(p-1))(1-u))^(q-1))du))
mean: (p/(p+q))
variance: ((p q)/((((p+q))^2)(p+q+1)))

r

dbeta(z,shape1 = p, shape2 = q)

# random number generation
rbeta(1,shape1 = p, shape2 = q)

mathematica

f[z_,p_,q_]:=(((z^(p-1))(((1-z))^(q-1)))/(Beta[p,q]))

Plot[f[z,2.2`,3.6`],{z,0,1}]

Therefore the support function is

(* beta distribution *)
sf[z_,p_,q_]:=Log[Beta[p,q]]-(p-1)Log[z]-(q-1)Log[(1-z)]

The Beta distribution can be transformed to be in terms of a mean parameter and variance parameter instead of p and q.  This can be more useful because one is often interested in modeling the mean.

substiution={q->((μ )/(σ^2))(μ(1-μ)-(σ^2))((1-μ )/μ),p->((μ )/(σ^2))(μ(1-μ)-(σ^2))};

f[z_,μ_,σ_]:=(((((1-z))^(-1+(((1-μ) ((1-μ) μ-(σ^2)))/(σ^2)))) (z^(-1+((μ ((1-μ) μ-(σ^2)))/(σ^2)))))/(Beta[((μ ((1-μ) μ-(σ^2)))/(σ^2)),(((1-μ) ((1-μ) μ-(σ^2)))/(σ^2))]));
sf[z_,μ_,σ_]:=(1-(((-1+μ) (-μ+(μ^2)+(σ^2)))/(σ^2)))Log[(1-z)]+(((-(μ^2)+(μ^3)+(σ^2)+μ (σ^2))/(σ^2)))Log[ z]+Log[Beta[-((μ (-μ+(μ^2)+(σ^2)))/(σ^2)),(((-1+μ) (-μ+(μ^2)+(σ^2)))/(σ^2))]] ;

where μ ∈ [0,1]
and (σ^2)<μ(1-μ)

mean: μ
variance: (σ^2)

The Beta distribution can be transformed to be in terms of a mean parameter and dispersion parameter, p = s(1-m) and q =s m.  This just scales p and q without the specific constraints on the variance version.  Here m ∈ (0,1), s ∈ (0,∞).

Bimodal Beta

extend the Beta Distribution by adding 2 pdf's together and dividing by 2
probability density function:
	f(z)=((((z^(p-1))(1-z))^(q-1))/(2B(p,q)))+((((z^(r-1))(1-z))^(s-1))/(2B(r,s))), B is the beta function, B(p,q) = (∫_0^1)(((u^(p-1))(1-u))^(q-1))du)
mean: 
variance:

mathematica

f[z_,v_,w_,x_,y_]:=(((z^(v-1))(((1-z))^(w-1)))/(2Beta[v,w]))+(((z^(x-1))(((1-z))^(y-1)))/(2Beta[x,y]))

f[0.4,2.2,3.6,2.2,3.6]

Plot[f[z,1.3`,3.6`,4.1`,1.6`],{z,0,1}]

Integrate[f[z,1.3,3.6,4.1,1.6],{z,0,1}]

Therefore the support function is

(* bimodal beta distribution *)
sf[z_,v_,w_,x_,y_]:=-Log[(((z^(v-1))(((1-z))^(w-1)))/(2Beta[v,w]))+(((z^(x-1))(((1-z))^(y-1)))/(2Beta[x,y]))]

Discrete

Expected, or mean
	E[Z]=(∑_z)z (Subscript[f, z])
Variance
	VAR(Z)=E[(((Z-E[Z]))^2)]=(∑_z)(((z-E[z]))^2)(Subscript[f, z])	
where Z is a random variable, z is a particular instance of that variable, (Subscript[f, z]) is the frequency of occurance of z.  

For discrete variables, the probability distribution function is the probability of the value z given the parameters

Binomial

Binomial Distribution:
has two parameters: N is the number of times the trial was conducted in which one of two outcomes can occur. p is the probability of one outcome, and z is the total number of times that outcome occurs.  For example, flip a coin 10 times and count the number of heads.  The result would be a binomial distribution, and if the coin is fair, p = 0. 5.  Whereas the other discrete distributions are used to model discrete observations of continuous processes, the binomial models discrete processes.
probability distribution function:
	Pr{Z=z}=(((N!)/((N-z)!z!)))(((p^z)(1-p))^(N-z))
mean: N p
variance: N p (1-p)

Note:  Sampling without replacment is best captured by a hypergeometric distribution

r

dbinom(z,n,p)
plot(dbinom(seq(1:20),20,0.3))

# random number generation
rbinom(1,n,p)

#see function description below
curve(1/(1+exp(-x)),-10,10)

mathematica

f[z_,n_,p_]:=((n!)/((n-z)!z!))(p^z)(((1-p))^(n-z))

ListPlot[Table[{z,f[z,20,0.25`]},{z,0,20}],PlotStyle->PointSize[0.03`]]

Sums to 1

Plus@@Table[f[z,20,0.3],{z,0,20}]

(* binomial random number with n = 3 and p = 0.5 *)
RandomInteger[BinomialDistribution[3 ,0.5]]

Here is a useful function for modeling a probability.  The funciton is constrained to a range between 0 and 1 for all real values of f(x), both positive and negative.  This function is derived by the following logic:

1) Probabilities ∈ [0,1].  An odds ratio is a ratio of the probability of sucess over the probability of failure, 
and odds ratios ∈ [0,∞).  Therefore
ln(odds ratios) ∈ (-∞,∞).

2)  For a binary model, only two outcomes are examined, therefore the Pr{failure} = 1 - Pr{success}.
ln(((Pr{success})/(Pr{failure})))=f(x)
((Pr{success})/(Pr{failure}))=exp(f(x))
((Pr{success})/(1-Pr{success}))=exp(f(x))
Pr{success}=((exp(f(x)))/(exp(f(x))+1))
Pr{success}=(1/(1+exp(-f(x))))

parameters={a->-4,b->0.7`};
f[x_]:=a+b x/.parameters;
Plot[(1/(1+Exp[-f[x]]))/.parameters,{x,0,10},PlotRange->{0,1}]

The Hill function is also useful for modeling probabilities as it contains interpretable parameters: a is the point at which the probability equals 0.50 and b is the "slope" at the 0.50 point.

parameters={a->5,b->4};
Plot[(1/(1+((((a/x)))^b)))/.parameters,{x,0,10},PlotRange->{0,1}]

Multinomial

Multinomial Distribution:
This is an extension of the binomial to m possible outcomes each with a probability (Subscript[p, i]).
probability distribution function:
	Pr{(Subscript[Z, m])=(Subscript[z, m])}=(((N!)/((Subscript[z, 1])!(Subscript[z, 2])! . . . (Subscript[z, m])!)))((Subscript[p, 1])^(Subscript[z, 1]))((Subscript[p, 2])^(Subscript[z, 2])) . . . ((Subscript[p, m])^(Subscript[z, m])) 
mean: N (Subscript[p, m])
variance: N (Subscript[p, m]) (1-(Subscript[p, m]))

mathematica

f[z1_,z2_,n_,p1_,p2_]:=((n!)/(z1!z2!(n-z1-z2)!))(p1^z1)(p2^z2)(((1-p1-p2))^(n-z1-z2))

ListPlot[Table[{z,f[z,3,20,0.25`,0.25`]},{z,0,20}],PlotStyle->PointSize[0.03`]]

Sums to 1

Plus@@Flatten[Table[f[z,i,5,0.25,0.25],{i,0,5},{z,0,5}]]

<<"MultivariateStatistics`"

(* multinomial random vector for an event that occurs 10 times; the probability vector follows *)
Random[MultinomialDistribution[10,{0.2,0.2,0.3,0.3}]]

r

dmultinom()

#random number generation

Poisson

Poisson Distribution:
has one parameter that can really be broken into two parameters: r is the rate constant and t is the time.
One description of a poisson distribution:  Imagine a process with a constant probability of occurance r.  The event can occur at any time, and repeatedly.  However, if we measure at discrete time intervals, and count the number of times the event has occurred for discrete entities, the result will be a Poisson Distribution.  Imagine rain falling onto floor tiles.  If you count the number of drops that hit each tile in one minute, the result would be distribute as a Poisson distribution.
probability distribution function:
	Pr{Z=z}=((((e^(-r t))(r t))^z)/(z!))
mean: r t
variance: r t
note: if evenly spaced time intervals, set t equal to 1

r

dpois(z,lambda)

# random number generation
rpois(1,lambda)

mathematica

f[z_,λ_]:=((Exp[-λ](λ^z))/(z!))

ListPlot[Table[{z,f[z,3]},{z,0,15}],PlotStyle->PointSize[0.03`]]

Sums to 1

Plus@@Table[f[z,20,0.3],{z,0,100}]

<<Statistics`DiscreteDistributions`
(* poisson random number with λ = 3.26 *)
Random[PoissonDistribution[3.26]]

Negative binomial

Negative Binomial Distribution:
The negative binomial is similar to the poisson except that the data are overdispersed.  As a consequence, the variance will always be larger than the mean.  The negative binomial can be derived assuming the rate constant in the poisson is no longer constant but is gamma distributed.  As k->∞, the distribution approaches a Poisson.  k is the shape parameter.
probability distribution function:
	Pr{Z=z}=((Γ(k+z))/(Γ(k)z!))((((k/(k+m))))^k)((((m/(k+m))))^z)
	where Γ is the gamma function, Γ(n)=(∫_0^∞)Exp(-t)(t^(n-1))dt
mean: m
variance: m+((m^2)/k)

Alternative form: used to describe the number the number of failures that occur before n sucesses have occured, where the probability of success is p.

r

# note, "mu =" is needed to specify this particular parameterization
dnbinom(z,size = k, mu= m)

rnbinom(1,size=k,mu=m)

mathematica

f[z_,k_,m_]:=((Gamma[k+z])/(Gamma[k]z!))((((k/(k+m))))^k)((((m/(k+m))))^z)

ListPlot[Table[{z,f[z,10,3.2`]},{z,0,15}],PlotStyle->PointSize[0.03`]]

ListPlot[Table[{z,f[z,0.26`,10]},{z,0,100}],PlotStyle->PointSize[0.03`]]

Sums to 1

Plus@@Table[f[z,10,3],{z,0,100}]//N

Derive expected value

Sum[z f[z,k,m],{z,0,∞}]

Derive variance

Sum[(((z-Sum[z f[z,k,m],{z,0,∞}]))^2) f[z,k,m],{z,0,∞}]//Simplify

(*
Mathematica 6 now contains a random number generator for the NegativeBinomial.  You have to use it in the following form so as to conform to the above./
 RandomInteger[NegativeBinomialDistribution[k ,(k/(k+m))]] *)

RandomInteger[NegativeBinomialDistribution[3 ,(3/(3+10))]]

Old versions of Mathematica did not contain a built-in random number generator for this form of the negative binomial.  While not very elegant, the function below can be used to randomly generate numbers from a negative binomial distribution.

NegBiRandom[m_,k_]:=Block[{funtable,x,i},funtable=Table[((Gamma[k+x] ((((k/(k+m))))^k) ((((m/(k+m))))^x))/(Gamma[k] x!)),{x,0,250}];cdf=Table[(∑_i^endsum)funtable[[i]],{endsum,1,Length[funtable]}];x=RandomReal[];i=1;While[x>cdf[[i]],i++;];i-1]

NegBiRandom[10,3]

(* negative binomial distribution *)
Clear[sf];
sf[z_,m_?NumberQ,k_?NumberQ]:=If[m>0&&k>0,-Log[Gamma[k+z]]+Log[z!]+Log[Gamma[k]]-k Log[k]+k Log[k+m]-z Log[m]+z Log[k+m],10000.0];
SetAttributes[sf,Listable];

Generalized Poisson

The generalized Poisson distribution (Consul and Jain. 1973. Technometrics. 15: 791-799) can be used to represent a Poisson process in which the data are over-dispersed.  It is also known as Lagrangian Poisson Distribution (Johnson, Kotz, and Kemp.  1992. Univariate Discrete Distributions).  Despite the claims of the original authors, this distribution cannot be used for under-dispersed data because the sum of the probabilities do not add to 1.

This distribution explicity models an increase in the rate of moving to subsequent categories and can represent data with very long right-tails as compared to the negative binomial.  This distribution worked very well for the distribution of pollen deposited on flowers (Castellanos et al. 2003.  Evolution 57:2742-2752.)  When (Subscript[λ, 2])=0 this distribution becomes the Poisson distribution.

Pr{Z=z}=(((e^(-((Subscript[λ, 1])+z (Subscript[λ, 2]))))(((Subscript[λ, 1])((Subscript[λ, 1])+z (Subscript[λ, 2])))^(z-1)))/(z!)), where  0<(Subscript[λ, 1]) and 0≤(Subscript[λ, 2])<1
mean: ((Subscript[λ, 1])/(1-(Subscript[λ, 2])))
variance: ((Subscript[λ, 1])/(((1-(Subscript[λ, 2])))^3))

f[z_,λ1_,λ2_]:=((λ1(((λ1+z λ2))^(z-1))Exp[-(λ1+z λ2)])/(z!))

Consul and Jain (1973) list the following method for estimating (Subscript[λ, 1])and (Subscript[λ, 2])from the mean and variance, but I don't trust these as the the best estimates, just a close estimate.

(Subscript[λ, 2])=1-(Sqrt[(mean/variance)])
(Subscript[λ, 1])=mean (1-(Subscript[λ, 2]))

Beta-Binomial

Beta Binomial Distribution:
This distribution assumes that p in the binomial distribution has a beta distribution with parameters v and w (see beta distribution above).

probability distribution function:
	Pr{Z=z}=(((N!)/((N-z)!z!)))((B(v+z,N+w-z))/(B(v,w))), where B(v,w) is the beta function
mean: N(p/(p+q))
variance: ((N p q (N+p+q))/((((p+q))^2)(1+p+q)))

The Beta distribution can be transformed to be in terms of a mean parameter and dispersion parameter, p = s(1-m) and q =s m.  This just scales p and q without the specific constraints on the variance version.  Here m ∈ (0,1), s ∈ (0,∞).  The m, s version is the one used for sfBetaBinomial in Mathematica.  This form has been used to describe correlated instances of the Binomial Distribution, although more mechanistic distributions for these cases exist.

mean: N m
variance: m (1-m)((N (N+s))/(1+s))
cv: ((√((N m(1 - m) (N + s))/(1 + s)))/(N m))

r

library(emdbook)
#this is based on  a parameterization so that m is the mean and s is 
#an overdispersion parameters.  As s approaches infintiy, 
#the Beta-Binomial collapses to a Biomial Distribution
dbetabinom(z,prob=m,size=N,theta=s)

mathematica

f[z_,n_,m_,s_]:=((n!)/((n-z)!z!))((Beta[s m+z,n+s (1-m)-z])/(Beta[s m,s (1- m)]))

f[z_,n_,v_,w_]:=((n!)/((n-z)!z!))((Beta[v+z,n+w-z])/(Beta[v,w]))

ListPlot[Table[{z,f[z,20,2.`,3.2`]},{z,0,20}],PlotStyle->PointSize[0.03`]]

Sums to 1

Plus@@Table[f[z,20,2.0,3.2],{z,0,20}]

(* Beta binomial distribution *)
sf[z_,n_,v_,w_]:=Log[(n-z)!]+Log[z!]-Log[n!]+Log[Beta[v,w]]-Log[Beta[v+z,n+w-z]]

The beta binomial is derived by multiply the pdf of the beta distribution by the pdf of the binomial, and then integrating from across the the full range of z for the beta binomial, in this case [0,1]

(* pdf for beta distribution *)
f[z_,p_,q_]:=(((z^(p-1))(((1-z))^(q-1)))/(Beta[p,q]))

(* pdf for binomial distribution *)
f[z_,n_,p_]:=((n!)/((n-z)!z!))(p^z)(((1-p))^(n-z))

Integrate[((n!)/((n-z)!z!))(p^z)(((1-p))^(n-z))(((p^(v-1))(((1-p))^(w-1)))/(Beta[v,w])),{p,0,1},Assumptions->{v>0,w>0,z≥0,n≥0}]

Mathematica doesn't do a good job of pulling out the beta function, but this simplifies to

((n!)/((n-z)!z!))((Beta[z+v]Beta[w+n-z])/(Beta[v,w]))

Bimodal Beta-Binomial

Using the same assumptions from the derivation of the unimodal beta binomial distribution, multiply the bimodal beta distribition by the binomial distribution and integrate p from 0 to 1

(* pdf for binomial distribution *)
f[z_,n_,p_]:=((n!)/((n-z)!z!))(p^z)(((1-p))^(n-z))

(* pdf for bimodal beta *)
f[z_,v_,w_,x_,y_]:=(((z^(v-1))(((1-z))^(w-1)))/(2Beta[v,w]))+(((z^(x-1))(((1-z))^(y-1)))/(2Beta[x,y]))

(* multiply both distributions by each other *)
((n!)/((n-z)!z!))(p^z)(((1-p))^(n-z))((((p^(v-1))(((1-p))^(w-1)))/(2Beta[v,w]))+(((p^(x-1))(((1-p))^(y-1)))/(2Beta[x,y])))
((n!)/((n-z)!z!))(p^z)(((1-p))^(n-z))(((p^(v-1))(((1-p))^(w-1)))/(2Beta[v,w]))+((n!)/((n-z)!z!))(p^z)(((1-p))^(n-z))(((p^(x-1))(((1-p))^(y-1)))/(2Beta[x,y]))
((n!)/((n-z)!z!))(((p^(z+v-1))(((1-p))^(w+n-z-1)))/(2Beta[v,w]))+((n!)/((n-z)!z!))(((p^(z+x-1))(((1-p))^(y+n-z-1)))/(2Beta[x,y]))
((n!)/((n-z)!z!))((((p^(z+v-1))(((1-p))^(w+n-z-1)))/(2Beta[v,w]))+(((p^(z+x-1))(((1-p))^(y+n-z-1)))/(2Beta[x,y])))
(∫_0^1)((n!)/((n-z)!z!))((((p^(z+v-1))(((1-p))^(w+n-z-1)))/(2Beta[v,w]))+(((p^(z+x-1))(((1-p))^(y+n-z-1)))/(2Beta[x,y])))\[DifferentialD]p
((n!)/((n-z)!z!))(∫_0^1)(((p^(z+v-1))(((1-p))^(w+n-z-1)))/(2Beta[v,w]))\[DifferentialD]p+((n!)/((n-z)!z!))(∫_0^1)(((p^(z+x-1))(((1-p))^(y+n-z-1)))/(2Beta[x,y]))\[DifferentialD]p
((n!)/((n-z)!z!))((Beta[z+v,w+n-z])/(2Beta[v,w]))+((n!)/((n-z)!z!))((Beta[z+x,y+n-z])/(2Beta[x,y]))
((n!)/((n-z)!z!))(((Beta[z+v,w+n-z])/(2Beta[v,w]))+((Beta[z+x,y+n-z])/(2Beta[x,y])))

(* pdf for bimodal beta-binomal distribution *)
f[z_,n_,v_,w_,x_,y_]:=((n!)/((n-z)!z!))(((Beta[z+v,w+n-z])/(2Beta[v,w]))+((Beta[z+x,y+n-z])/(2Beta[x,y])))

Bimodal Beta-Binomial Distribution:
has 5 parameters: shape parameters for one node are v and w; shape parameters for the other node are x and y;  n is the total number of individuals in the set and z is the number of success in that set.

Note:  this distribution and the unimodal beta-binomial assume a different draw from beta for the p of each set, but it assumes the same p for all the individuals within a set.  This statement is based on Skellam. 1948. Journal of the Royal Statistical Society.  Series B, 10: 257-261.

probability density function:
	f(z)=((((z^(p-1))(1-z))^(q-1))/(B(p,q))), B is the beta function, B(p,q) = (∫_0^1)(((u^(p-1))(1-u))^(q-1))du)
	f(z)=((((z^(p-1))(1-z))^(q-1))/((∫_0^1)(((u^(p-1))(1-u))^(q-1))du))
mean: 
variance:

ListPlot[Table[f[z,30,1.2`,3.6`,3.6`,1.2`],{z,0,30}],PlotStyle->{PointSize[0.02`]},PlotRange->{0,0.04`}]

Plot[(((z^(v-1)) (((1-z))^(w-1)))/(2 Beta[v,w]))+(((z^(x-1)) (((1-z))^(y-1)))/(2 Beta[x,y]))/.{v->1.2`,w->3.6`,x->3.6`,y->1.2`},{z,0,1}]

Plus@@Table[f[z,30,1.2,3.6,3.6,1.2],{z,0,30}]

Therefore the support function is

(* bimodal beta-binomial distribution *)
sf[z_,n_,v_,w_,x_,y_]:=Log[(n-z)!]+Log[z!]-Log[n!]-Log[((Beta[z+v,w+n-z])/(2Beta[v,w]))+((Beta[z+x,y+n-z])/(2Beta[x,y]))]

Bimodal Binomial

f[z_,n_,p_,q_]:=((n!)/(2(n-z)!z!))(p^z)(((1-p))^(n-z))+((n!)/(2(n-z)!z!))(q^z)(((1-q))^(n-z))

ListPlot[Table[{z,f[z,20,0.15`,0.45`]},{z,0,20}],PlotStyle->PointSize[0.03`]]

Sums to 1

Plus@@Table[f[z,20,0.3],{z,0,20}]

(* bimodal binomial *)
sf[z_,n_,p_,q_]:=Log[2(n-z)!]+Log[z!]-Log[n!]-Log[(p^z)(((1-p))^(n-z))+(q^z)(((1-q))^(n-z))]

Load data

This section will load data from a "csv" (Comma delimited) file in which each row is a replicate.  Microsoft Excel can save worksheets in "csv" format.  Each column should have a text heading.
	-start each heading with a small letter
	-do not have any spaces or periods in the column heading
	-the values for categorical variables need to be strings, not numbers
	
Modify the location in the following line.

data=Import["C:\\Users\\Chad\\Documents\\webs\\Chad Brassil\\ELME\\2007\\eejData.csv","CSV"];

Alternatively. you can manually enter data.  You can insert a matrix template by right-clicking and choosing "Create Table/Matrix/Palette..."  The top row should be column headings with the text in quotation marks.  Categorical variables should also be in quotation marks.

(* example of manual data entry *)
(*data=(({{"size", "frogs", "species"}, {0.23, 1, "a"}, {0.51, 0, "b"}}));*)

The ProcessData[ ] command will attempt to extract column headings from the data file.  The column headings will be placed in "labels" and the remaining numbers will be left in "data".  Categorical variables will be identified and dummy variables (i.e. 1 or 0) will be generated for each category.

ProcessData[data]

The following lines display the data and the column labels

MatrixForm[data]

labels

special functions

The CPlot[ ] function can be used to graph data from two columns

(* enter column labels *)
pObs=CPlot[mInfl,visits]

The SelectCases[ ] function can be used to select a subset of the data based on a categorical variable

SelectCases[mInfl,species,a]

Use list to select multiple categories

SelectCases[mInfl,species,{a,b}]

CPlot3D plots a 3D scatter plot

CPlot3D[mInfl,knInfl,visits];

General maximum likelihood code

In addition, to these steps, one should first build a model and simulate the data

(* First, load FinalTable where results will be stored *)
FinalTable={};

Choose distribution for support function

Choose your support function based on the probability density function, pre-build support functions are listed here

?sfPoisson
?sfNegBinom
?sfGenPoisson
?sfBinomial
?sfNormal
?sfBimodalNormal
?sfLogNormal
?sfBetaMean
?sfBetaMeanShape
?sfBetaBinomial

Model A: maximum likelihood estimates

Customize the negative-log-likelihood function (ll) to reflect your model and your support function

Clear[ll];
ll[]:=Plus@@
sfNormal[visits,β0+β1 mInfl,σ];

You should be able to use the following code without modification.  If the likelihood surface is particularily complex, try alternative minization functions below.  Constraints can be manually entered, see alternative minization functions below.

vars=ProcessVars[ll[]];constraints=ProcessConstraints[ll];
ans=NMinimize[{ll[],constraints},vars][[2]]

maxll=ll[]/.ans

modelName="constant";
AppendTo[FinalTable,{maxll,Length[vars],modelName}];

Plot model and data

pPred=Plot[β0+β1 mInfl/.ans,{mInfl,Min[mInfl],Max[mInfl]}]
Show[pObs,pPred]

Categorical Examples

Use the function Categorical to create seperate parameters for each category in a categorical variable.  In this case, individual β0 parameters will be fit for each species.

Clear[ll];
ll[]:=Plus@@
sfNormal[visits,Categorical[β0,species]+β1 mInfl,σ];

vars=ProcessVars[ll[]];constraints=ProcessConstraints[ll];
ans=NMinimize[{ll[],constraints},vars][[2]]

pObs=CPlot[SelectCases[mInfl,species,a],SelectCases[visits,species,a],DisplayFunction->Identity];
pPred=Plot[(Subscript[β0, a])+β1 mInfl/.ans,{mInfl,Min[mInfl],Max[mInfl]},DisplayFunction->Identity];
Show[pObs,pPred,DisplayFunction->$DisplayFunction]

pObs=CPlot[SelectCases[mInfl,species,b],SelectCases[visits,species,b],DisplayFunction->Identity];
pPred=Plot[(Subscript[β0, b])+β1 mInfl/.ans,{mInfl,Min[mInfl],Max[mInfl]},DisplayFunction->Identity];
Show[pObs,pPred,DisplayFunction->$DisplayFunction]

Categories can also be combined.  In this case, species groups a and b are combined.

Clear[ll];
ll[]:=Plus@@
sfNormal[visits,Categorical[β0,species,{a->ab,b->ab}]+β1 mInfl,σ];

vars=ProcessVars[ll[]];constraints=ProcessConstraints[ll];
ans=NMinimize[{ll[],constraints},vars][[2]]

Alternative  Functions (if needed)

(* can specify parameter constraints directly *)
constraints=a>0&&b>0;
ans=NMinimize[{ll[],constraints},vars][[2]]

(* can specify inital rectangular regions *)
upperguess={12.0,8.0};
lowerguess={7.0,3.0};
ans=NMinimize[{ll[],constraints},Transpose[{vars,lowerguess,upperguess}]][[2]]

(* code that will work with NMinimize[], internally this uses FindMinimum from random starting positions *)
NMinimize[{ll[],constraints},vars,Method->"RandomSearch"][[2]]

guess={50.0,7.0,3.0}; (* initial guess, not used with NMinimize, but used with FindRoot and FindMinimum *)

(* standard form of FindRoot uses Newton methods, but now NMinimize is better because it uses Nelder-Mead methods *)
derv=Map[D[ll[],#]&,vars];
ans=FindRoot@@Join[{Thread[derv==0]},Transpose[{vars,guess}]]

(* slightly different form of FindRoot that you can also try *)
ans=FindRoot@@Join[{Thread[derv==0]},Transpose[{vars,guess,guess*0.98}]]

(* if FindRoot[] does not work, can also try FindMinimum[] *)
ans=FindMinimum[ll[],Transpose[{vars,guess}]][[2]]

guessRange={{0.`,3.`},{0.`,2.`},{0,1},{0,0.1`}};
seeds=Table[RandomReal[guessRange[[i]]],{1000},{i,Length[guessRange]}];
seedRes=Apply[ll,seeds,1];
guess=seeds[[Position[seedRes,Min[Cases[seedRes,_Real]]][[1,1]]]]
ll@@guess

(* more explicit ll function that might be needed for specialized data *)
ll[m_,b_,σ_]:=Plus@@
Table[
sf[data[[i]],m+b*i,σ]
,{i,Length[data]}];

(* can specify vars list directly *)
vars={a,b,σ};

Confidence Intervals

This function is for use with the NMinimize set-up.

This code will automatically attempt to find confidence intervals for all functions.  Note, you first need to execute llProfileFun[] above.  It may take a long time, especially if it has problems with particular parameters.

llProfileFun[paramNo_,paramValue_]:=Module[{tempVars,newans},
tempVars=Drop[vars,{paramNo}];
newans=NMinimize[{ll[],constraints}/.vars[[paramNo]]->paramValue,tempVars][[2]];
ReplaceAll[ll[],Join[newans,{vars[[paramNo]]->paramValue}]]
];
ci={};
ci=Do[
AppendTo[ci,{paraMin,paraMax}/.FindRoot[Hold[llProfileFun[varNo,paraMin]==maxll+1.92],{paraMin,vars[[varNo]]/.ans,(vars[[varNo]]/.ans)*0.9}]/.FindRoot[Hold[llProfileFun[varNo,paraMax]==maxll+1.92],{paraMax,vars[[varNo]]/.ans,(vars[[varNo]]/.ans)*1.1}]]
,{varNo,Length[vars]}];
MatrixForm[Transpose[{vars,vars/.ans,ci}]]

(* check to see if c.i. can be reported as mean ± c.i. *)
MatrixForm[Transpose[{vars,Table[{ci[[i,1]]-ans[[i,2]],ci[[i,2]]-ans[[i,2]]},{i,Length[vars]}]}]]

For use with alternative functions

(* If using lowerguess and upperguess vectors to constrain NMinimize *)

(* find confidence intervals using likelihood profile *)
llProfileFun[paramNo_,paramValue_]:=Module[{tempVars,newans},
tempVars=Drop[vars,{paramNo}];
tempUpperguess=Drop[upperguess,{paramNo}];
tempLowerguess=Drop[lowerguess,{paramNo}];
newans=NMinimize[{ll@@vars}/.vars[[paramNo]]->paramValue,Transpose[{tempVars,tempLowerguess,tempUpperguess}]][[2]];
Apply[ll,ReplaceAll[vars,Join[newans,{vars[[paramNo]]->paramValue}]]]
];
ci=Table[
{paraMin,paraMax}/.FindRoot[Hold[llProfileFun[varNo,paraMin]==maxll+1.92],{paraMin,vars[[varNo]]/.ans,(vars[[varNo]]/.ans)*0.9}]/.FindRoot[Hold[llProfileFun[varNo,paraMax]==maxll+1.92],{paraMax,vars[[varNo]]/.ans,(vars[[varNo]]/.ans)*1.1}]
,{varNo,Length[vars]}];
MatrixForm[Transpose[{vars,vars/.ans,ci}]]

(* For use with alt. min. functions *)

(* find confidence intervals using likelihood profile *)
llProfileFun[paramNo_,paramValue_]:=Module[{tempVars,tempGuess,derv,newans},
tempVars=Drop[vars,{paramNo}];
tempGuess=Drop[vars/.ans,{paramNo}];
derv=Map[D[Apply[ll,vars/.vars[[paramNo]]->paramValue],#]&,tempVars];
newans=FindRoot@@Join[{Thread[derv==0]},Transpose[{tempVars,tempGuess}]];
Apply[ll,ReplaceAll[vars,Join[newans,{vars[[paramNo]]->paramValue}]]]
];
ci=Table[
{paraMin,paraMax}/.FindRoot[Hold[llProfileFun[varNo,paraMin]==maxll+1.92],{paraMin,vars[[varNo]]/.ans,(vars[[varNo]]/.ans)*0.9}]/.FindRoot[Hold[llProfileFun[varNo,paraMax]==maxll+1.92],{paraMax,vars[[varNo]]/.ans,(vars[[varNo]]/.ans)*1.1}]
,{varNo,Length[vars]}]

Optional: examine likelihood surface for a parameter

problemParameter=1;
problemRange={3.5`,4.5`};
numberOfDivisions=15;

llProfileFun[paramNo_,paramValue_]:=Module[{tempVars,newans},tempVars=Drop[vars,{paramNo}];newans=NMinimize[{ll[],constraints}/.vars[[paramNo]]->paramValue,tempVars][[2]];ll[]/.Join[newans,{vars[[paramNo]]->paramValue}]];
SetAttributes[llProfileFun,Listable];
sampleAt=Table[i,{i,problemRange[[1]],problemRange[[2]],((problemRange[[2]]-problemRange[[1]])/(numberOfDivisions-1))}];

Show[ListPlot[Transpose[{sampleAt,llProfileFun[problemParameter,sampleAt]}],PlotStyle->PointSize[0.02`],DisplayFunction->Identity],Plot[maxll+1.92`,{i,problemRange[[1]],problemRange[[2]]},DisplayFunction->Identity],DisplayFunction->$DisplayFunction]

Contour Plots of Confidence Intervals (using likelihood profile)

contourFunction[r_,k_]:=OneSidedPValue/.ChiSquarePValue[2 (ll[r,k]-maxll),2];
ContourPlot[contourFunction[r,k],{k,100,10000},{r,0.09`,0.17`},PlotPoints->50,Contours->{0.05`}]

OLD: Confidence Intervals with dragging starting conditions for problem parameters (using likelihood profile)

Clear[tempGuess];
llProfileFunDragGuess[paramNo_,paramValue_]:=Module[{tempVars,derv,newans},
tempVars=Drop[vars,{paramNo}];
derv=Map[D[Apply[ll,vars/.vars[[paramNo]]->paramValue],#]&,tempVars];
newans=FindRoot@@Join[{Thread[derv==0]},Transpose[{tempVars,tempGuess}]];
tempGuess=tempVars/.newans;
Apply[ll,ReplaceAll[vars,Join[newans,{vars[[paramNo]]->paramValue}]]]
];

varNo=1
tempGuess=Drop[vars/.ans,{varNo}];
FindRoot[Hold[llProfileFun[varNo,paraMin]==maxll+1.92],{paraMin,vars[[varNo]]/.ans,(vars[[varNo]]/.ans)*0.9}]
tempGuess=Drop[vars/.ans,{varNo}];
FindRoot[Hold[llProfileFun[varNo,paraMax]==maxll+1.92],{paraMax,vars[[varNo]]/.ans,(vars[[varNo]]/.ans)*1.1}]

OLD: Automatic Confidence Intervals with dragging starting conditions for problem parameters

Clear[tempGuess];
llProfileFunDragGuess[paramNo_,paramValue_]:=Module[{tempVars,derv,newans},
tempVars=Drop[vars,{paramNo}];
derv=Map[D[Apply[ll,vars/.vars[[paramNo]]->paramValue],#]&,tempVars];
newans=FindRoot@@Join[{Thread[derv==0]},Transpose[{tempVars,tempGuess}]];
tempGuess=tempVars/.newans;
Apply[ll,ReplaceAll[vars,Join[newans,{vars[[paramNo]]->paramValue}]]]
];

stepVector=N[10^(Floor[Log[10,Abs[(vars/.ans)]]]-1)]

ciResults={};
Do[
tempGuess=Drop[vars/.ans,{varNo}];
currentMaxll=maxll;
val=Take[vars/.ans,{varNo}];
While[currentMaxll<maxll+1.92,
previousVal=val;
previousMaxll=currentMaxll;
val+=stepVector[[varNo]];
currentMaxll=llProfileFunDragGuess[varNo,val];
];
middleVal=Round[(val/(stepVector[[varNo]]))]*stepVector[[varNo]]-((stepVector[[varNo]])/2);
middleMaxll=llProfileFunDragGuess[varNo,middleVal];
AppendTo[ciResults,{varNo,{previousVal,previousMaxll},{middleVal,middleMaxll},{val,currentMaxll}}];

tempGuess=Drop[vars/.ans,{varNo}];
currentMaxll=maxll;
val=Take[vars/.ans,{varNo}];
While[currentMaxll<maxll+1.92,
previousVal=val;
previousMaxll=currentMaxll;
val-=stepVector[[varNo]];
currentMaxll=llProfileFunDragGuess[varNo,val];
];
middleVal=Round[(val/(stepVector[[varNo]]))]*stepVector[[varNo]]+((stepVector[[varNo]])/2);
middleMaxll=llProfileFunDragGuess[varNo,middleVal];
AppendTo[ciResults,{varNo,{previousVal,previousMaxll},{middleVal,middleMaxll},{val,currentMaxll}}];,
{varNo,Length[vars]}]

maxll+1.92

MatrixForm[ciResults]

Table[{vars[[i]],Round[((ans[[i,2]])/(stepVector[[i]]))]*stepVector[[i]],If[ciResults[[2*i,3,2]]>maxll+1.92,Round[((ciResults[[2*i,2,1,1]])/(stepVector[[i]]))]*stepVector[[i]],Round[((ciResults[[2*i,4,1,1]])/(stepVector[[i]]))]*stepVector[[i]]],If[ciResults[[2*i-1,3,2]]>maxll+1.92,Round[((ciResults[[2*i-1,2,1,1]])/(stepVector[[i]]))]*stepVector[[i]],Round[((ciResults[[2*i-1,4,1,1]])/(stepVector[[i]]))]*stepVector[[i]]]},{i,Length[vars]}]//MatrixForm

OLD: Alternative Minimization Function (using likelihood surface)

problemParameter=1;
problemRange={7.`,8.`};
numberOfDivisions=10;

largestLL=2;
smallestLL=1;

While[largestLL-smallestLL>0.001`,sampleAt=Table[i,{i,problemRange[[1]],problemRange[[2]],((problemRange[[2]]-problemRange[[1]])/numberOfDivisions)}];tempVars=Drop[vars,{problemParameter}];tempGuess=Drop[guess,{problemParameter}];likelihoodProfile={};bestfits={};Do[derv=((∂_(#1))ll@@(vars/.vars[[problemParameter]]->sampleAt[[i]])&)/@tempVars;ans=FindRoot@@Join[{Thread[derv==0]},Transpose[{tempVars,tempGuess}]];AppendTo[bestfits,ans];AppendTo[likelihoodProfile,{sampleAt[[i]],ll@@vars/.vars[[problemParameter]]->sampleAt[[i]]/.ans}];,{i,numberOfDivisions}];Print[ListPlot[likelihoodProfile,PlotStyle->PointSize[0.02`]]];largestLL=Max[Column[likelihoodProfile,2]];smallestLL=Min[Column[likelihoodProfile,2]];pos=Position[Column[likelihoodProfile,2],smallestLL][[1,1]];guess=vars/.bestfits[[pos]]/.vars[[problemParameter]]->likelihoodProfile[[pos,1]];problemRange={likelihoodProfile[[pos,1]]-((problemRange[[2]]-problemRange[[1]])/numberOfDivisions),likelihoodProfile[[pos,1]]+((problemRange[[2]]-problemRange[[1]])/numberOfDivisions)};];
smallestLL
Thread[vars->guess]

Akaike's Information Criterion (AIC)

Compares nested or non-nested models.  The total number of parameters is n.  The model with the lowest AIC is the most likely model.
AIC = 2 maxll + 2 parameters

AICtable=Transpose[Append[Transpose[FinalTable],2 FinalTable[[All,1]]+2 FinalTable[[All,2]]]];

minAIC=Min[AICtable[[All,4]]];
AICtable=Transpose[Append[Transpose[AICtable],AICtable[[All,4]]-minAIC]];


AICtable=Transpose[Append[Transpose[AICtable],((Exp[-0.5` AICtable[[All,5]]])/(Plus@@Exp[-0.5` AICtable[[All,5]]]))]];


AICtable=Transpose[Append[Transpose[AICtable],Exp[0.5` AICtable[[All,5]]]]];
MatrixForm[Prepend[Sort[AICtable,#1[[7]]<#2[[7]]&],(Style[#1,FontWeight->"Bold"]&)/@{ll,parameters,model,AIC,(Subscript[Δ, i]),Akaike weights,evidence ratio}]]

corrected Akaike's Information Criterion ((Subscript[AIC, c]))

The (Subscript[AIC, c]) values correct for small sample sizes because AIC values perform poorly when a large number of parameters are being estimated relative to the sample size.  Generally, the (Subscript[AIC, c]) is highly recommended when the ratio of ((sample size)/parameters)< 40  (Burnham and Anderson 2002).

(*** MAY NEED TO MODIFY THIS LINE ***)
sampleSize=Length[data]
(*** MAY NEED TO MODIFY THIS LINE ***)

AICtable=Transpose[Append[Transpose[FinalTable],2 FinalTable[[All,1]]+2 FinalTable[[All,2]]]];

AICtable=Transpose[Append[Transpose[AICtable],AICtable[[All,4]]+((2 FinalTable[[All,2]] (FinalTable[[All,2]]+1))/(sampleSize-FinalTable[[All,2]]-1))]];

minAICc=Min[AICtable[[All,5]]];
AICtable=Transpose[Append[Transpose[AICtable],AICtable[[All,5]]-minAICc]];


AICtable=Transpose[Append[Transpose[AICtable],((Exp[-0.5` AICtable[[All,6]]])/(Plus@@Exp[-0.5` AICtable[[All,6]]]))]];


AICtable=Transpose[Append[Transpose[AICtable],Exp[0.5` AICtable[[All,6]]]]];
MatrixForm[Prepend[Sort[AICtable,#1[[8]]<#2[[8]]&],(Style[#1,FontWeight->"Bold"]&)/@{ll,parameters,model,AIC,AICc,(Subscript[Δ, i]),Akaike weights,evidence ratio}]]

Likelihood ratio test: for nested models

Compares two nested models having likelihoods maxllA nad maxllB.  Degrees of freedom, df, is the number of parameters that differ between the models.  The resulting p-value is interpretted as in classical statistics.  If p < 0.05, we reject the hypothesis that there is no difference between the models.  In other words, if p < 0.05, we can say the more complex model is a better fit.

Table[Append[FinalTable[[i]],i],{i,Length[FinalTable]}]//MatrixForm

compA=5; (* model number in Final Table *)
compB=4;(* model number in Final Table *)

df=Abs[FinalTable[[compA,2]]-FinalTable[[compB,2]]];
r=Abs[2*(FinalTable[[compA,1]]-FinalTable[[compB,1]])];
ChiSquarePValue[r,df]

Manual Comparison

df=1;
r=2(maxllA-maxllB)
ChiSquarePValue[r,df]