{ From user Lonnie, Model Multivariate_distrib at Tue, Nov 11, 2008 1:59 PM~~
}
Softwareversion 4.2.0



Linklibrary Multivariate_distrib
Title: Multivariate Distributions
Description: A library of multivariate distributions.~
~
In a multivariate distribution, each sample is a vector.  This vector~~
 is identified by an index, identified by the I parameter of the func~~
tions in this library.  A Mid value from a distribution function will~~
 therefore be indexed by I, whlie a Sample from a distribution functi~~
on is indexed by both I and Run.  These distribution functions can al~~
so be used from within the Random function to generate a single monte~~
-carlo sample, which will be indexed by I.~
~
This library also contains functions for generating correlated distri~~
butions.  Correlate_with, for example, allows you to generate a univa~~
rite distribution with an arbitrary marginal distribution that has a ~~
specified rank correlation with an arbitrary reference distribution. ~~
  Several functions may be used for generating serial correlations, w~~
here each distribution along an index is correlated with the previous~~
 point along that index.
Author: Lonnie Chrisman, Ph.D.~
Lumina Decision Systems~
~
With contributions by:~
   John Bowers, US FDA.~
   Max Henrion, Lumina Decision Systems
Date: Fri, Aug 01, 2003 7:12 PM
Saveauthor: Lonnie
Savedate: Tue, Nov 11, 2008 1:59 PM
Defaultsize: 48,24
Nodesize: 56,24
Nodeinfo: 1,1,1,1,1,1,0,0,0,0
Diagstate: 1,42,10,649,1009,17
Windstate: 2,401,199,483,316
Fontstyle: Arial, 15
Fileinfo: 0,Linklibrary Multivariate_distrib,2,2,0,0,W:\Analytica\Exec~~
Debug\Libraries\Multivariate Distributions.ana

Function Wishart( cv : Number[I,J,Run] ; n :positive ; I,J : Index ; ~
singleSampleMethod : optional hidden scalar)
Title: Wishart(cv,n,I,J)
Description: Suppose you sample N samples from a Gaussian(0,cv,I,J) di~~
stribution, X[I,R].  (R is the index that indexes each sample, R:=1..~~
N).  The Wishart distribution describes the distribution of sum( X * ~~
X[I=J], R ).  This matrix is dimensioned by I and J and is called the~~
 scatter matrix. ~
~
A sample drawn from the Wishart is therefore a sample scatter matrix.~~
  If you divide that sample by (N-1), you have a sampled covariance m~~
atrix.  ~
~
If you compute a sample covariance matrix from data, and then want to~~
 use this in your model, if you just use it directly, you'll be ignor~~
ing sampling error.  That may be insignificant of N is large.  Otherw~~
ise, you may want to use:~
  Wishart( SampleCV, N, I, J) / (N-1)~
instead of just SampleCV in your model.   The extended variance will ~~
account for the uncertainty from the finite sample size that was used~~
 to obtain your sample CV.~
~
If you can express a prior probability on covariances in the form of ~~
an InvertedWishart distribution, then the posterior distribution, aft~~
er having computed the sample covariance matrix (assumed to be drawn,~~
 by nature, from a Wishart), is also an InvertedWishart.
Definition: var T := if i<j then 0 else if i=j ~
             then sqrt(ChiSquared(n-i+1, singleSampleMethod:singleSam~~
pleMethod))~
             else (normal(0,1, singleSampleMethod:singleSampleMethod)~~
);~
var A := MatrixMultiply(T,I,J,Transpose(T,I,J),I,J);~
var L := Decompose(A,I,J);~
MatrixMultiply(MatrixMultiply(L,I,J,A,I,J),I,J,Transpose(L,I,J),I,J)~~

Nodelocation: 128,224,1
Nodesize: 100,16
Windstate: 2,183,11,538,656
Paramnames: cv,n,I,J

Function Invertedwishart(psi:Number[I,J,Run] ; n :positive ; I,J : Ind~~
ex ;~
singleSampleMethod : optional hidden scalar)
Title: Inverted Wishart(Psi,n,I,J)
Description: The inverted Wishart distribution represents a distributi~~
on of covariance matrices.  Each sample from the InvertedWishart is a~~
 covariance matrix.  It is conjugate to a Wishart distribution, and t~~
herefore is sometimes used as a Bayesian prior distribution for covar~~
iance.  The parameter, Psi, must be a positive definite matrix.~
~
Suppose you represent the prior distribution of covariance using an i~~
nverted Wishart distribution: InvertedWishart(Psi,m).  You observe so~~
me data, X[I,R], where R:=1..N indexes each datapoint and I is the ve~~
ctor dimension, and compute A = Sum( X*X[I=J], R), where A is called ~~
the scatter matrix.  A is an observation that gives you information a~~
bout the true covariance matrix, so can use this to obtain a Bayesian~~
 posterior distribution on the true covariance given by:~
  InverseWishart( A+Psi, n+m ).~

Definition: Invert( Wishart( Invert(Psi,I,J), n, I,J, singleSampleMeth~~
od:singleSampleMethod), I, J )
Nodelocation: 352,224,1
Nodesize: 108,16
Windstate: 2,102,90,553,530
Paramnames: psi,n,I,J

Function Regressiondist( Y : Numeric[I,Run] ; B : Numeric[I,K,Run] ; I~~
,K : Index; C : optional Numeric[K,Run] ;~
 S : optional Numeric[I,Run] ;~
singleSampleMethod : optional hidden scalar )
Title: RegressionDist(Y,B,I,K)
Description: RegressionDist returns linear regression coefficients as ~~
a distribution.~
~
Suppose you have data where Y was produced as:~
  Y = Sum( C*B, K ) + Normal(0,S)~
~
S is the measurement noise.  You have the data (B[I,K] and Y[I]).  Yo~~
u might or might not know the measurement noise S.  So you perform a ~~
linear regression to obtain an estimate of C.  Because your estimate ~~
is obtained from a finite amount of data, your estimate of C is itsel~~
f uncertain.  This function returns the coefficients C as a distribut~~
ion (i.e., in Sample mode, it returns a sampling of coefficients inde~~
xed by Run and K), reflecting the uncertainty in the estimation of th~~
ese parameters.~
~
If you know the noise level S in advance, then you can use historical~~
 data as a starting point for building a predictive model of Y, as fo~~
llows:~
~
{ Your model of the dependent variables: }~
Variable Y := your historical dependent data, indexed by I~
Variable B := your historical independent data, indexed by I,K~
Variable X := { indexed by K.  Maybe others.  Possibly uncertain }~
Variable S := { the known noise level }~
Chance C := RegressionDist(Y,B,I,K)~
Variable Predicted_Y := Sum(C*X,K) + Normal(0,S)~
~
If you don't know the noise level, then you need to estimate it. You'~~
ll need it for the normal term of Predicted_Y anyway, and you'll need~~
 to do a regression to find it.  So you can pass these optional param~~
eters into RegressionDist.  The last three lines above become:~
Variable E_C := Regression(Y,B,I,K)~
Variable S := RegressionNoise( Y,B,I,K,E_C )~
Chance C := RegressionDist(Y,B,I,K,E_C)~
Variable Predicted_Y := Sum(C*X,K) + Normal(0,S)~
~
If you use RegressionNoise to compute S, you should use Mid(Regressio~~
nNoise(...)) for the S parameter.  However, when computing S for your~~
 prediction, don't RegressionNoise in context.  Better is if you don'~~
t know the measurement noise in advance, don't supply it as a paramet~~
er.
Definition: if IsNotSpecified(C) Then C := Regression(Y,B,I,K);~
If IsNotSpecified(S) Then S := Mid(RegressionNoise(Y,B,I,K,C))~
                                 Else S:=Mean(S);~
Index K2 := K;~
var cv := Invert( Sum( B*B[K=K2]/S^2,I),K,K2);~
cv := Average( [cv,Transpose(cv,K,K2)]); {for round off errs}~
~
Gaussian(C,cv,K,K2,singleSampleMethod:singleSampleMethod)
Nodelocation: 152,736,1
Nodesize: 112,20
Windstate: 2,115,1,561,882
Paramnames: Y,B,I,K,C,S

Function Gaussian(meanVec : numeric[I,Run],covar : numeric[I,J,Run]; I~~
,J:IndexType ;~
Over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Gaussian( m, cv, i, j )
Description: A multi-variate Gaussian distribution based on a mean vec~~
tor and covariance matrix.  The covariance matrix must symmetric and ~~
positive-definite.  The meanVec is indexed by I.  The covariance matr~~
ix is 2-D, indexed by I & J.  Indexes I & J should be the same length~~
.
Definition: var S := Decompose(covar,I,J);~
var Z := Normal(0,1,over:J,singleSampleMethod:singleSampleMethod);~
sum( S*Z,J ) + meanVec
Nodelocation: 112,72,1
Nodesize: 80,16
Windstate: 2,96,299,486,314
Paramnames: meanVec,covar,I,J,Over

Function Dirichlet(alpha : Numeric[I,Run]; I:IndexType ;~
Over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar)
Title: Dirichlet ( a, i )
Description: A Dirichlet distribution with parameters alpha_i>0~
Each sample of a Dirichlet distribution produces a random vector whos~~
e elements sum to 1.  It is commonly used to represent second order p~~
robability information.~
~
The Dirichlet distribution has a density given by ~
   k * Product( X^(alpha-1), I)~
where k is a normalization factor equal to~
  GammaFn( sum(alpha,I )) / Sum(GammaFn(alpha),I)~
~
The parameters, alpha, can be interpreted as observation counts.  The~~
 mean is given by the relative values of alpha (normalized to 1), but~~
 the variance narrows as the alphas get larger, just as your confiden~~
ce in a distribution would narrow as you get more samples.~
~
The Dirichlet lends itself to easy Bayesian updating.  If you have a ~~
prior of alpha0, and you observe N
Definition: var a:=Gamma(alpha,singleSampleMethod:singleSampleMethod);~~
~
a/sum(a,I)
Nodelocation: 272,120,1
Nodesize: 58,16
Windstate: 2,26,18,624,485
Paramnames: alpha,I,Over

Function Binormal(MeanVec :numeric[I,Run]; Sdeviations : positive[I,Ru~~
n]; I:IndexType; correlationCoef : numeric[Run];~
Over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar)
Title: BiNormal (m, s, i, c )
Description: A 2-D Normal (or Bi-variate Gaussian) distribution with t~~
he indicated individual standard deviations (>0) and the indicated co~~
rrelation coefficient.  The index, I, must have exactly 2 elements, S~~
deviations must be indexed by I.
Definition: if size(I)<>2 then ~
   Error("Index to BiNormal must have 2 elements")~
else begin~
   var s := product(Sdeviations,I) * correlationCoef;~
   Index J:=CopyIndex(I);~
   Gaussian( meanVec, If I<>J Then s else Sdeviations^2, I,J,~
                        singleSampleMethod: singleSampleMethod )~
end
Nodelocation: 288,72,1
Nodesize: 78,16
Windstate: 2,2,24,525,540
Paramnames: MeanVec,Sdeviations,I,correlationCoef,Over

Function Multinomial(N:NonNegative ; theta:NonNegative  ; I : IndexTyp~~
e;~
Over : ... optional atomic ;~
singleSampleMethod : hidden optional scalar )
Title: Multinomial (n, theta, i )
Description: Returns the Multinomial Distribution.~
~
The multinomial distribution is a generalization of the Binomial dist~~
ribution to N possible outcomes.  For example, if you were to roll a ~~
fair die N times, an outcome would be the number of times each of the~~
 six numbers appears.  Theta would be the probability of each outcome~~
, where sum(theta,I)=1, and index I is the list of possible outcome. ~~
 If theta doesn't sum to 1, it is normalized.~
~
Each sample is a vector indexed by I indicating the number of times t~~
he corresponding outcome (die number) occurred during that sample poi~~
nt.  Each sample will have the property that sum( result, I ) = N.
Definition: var z := n;~
var k := size(I);~
~
var j:=cumulate(1,I) in I do begin~
  Index I2 := j..k;~
  var theta2 := Slice(theta,I,I2); /* unnormalized sub-process */~
  var p := theta2/sum(theta2, I2);~
  p := if IsNan(p) then 0 else p;~
  var xj := Binomial(z,p[I2=j],~
                              singleSampleMethod:singleSampleMethod);~~
~
  z := z - xj;~
  xj~
end~

Nodelocation: 117,120,1
Nodesize: 85,16
Windstate: 2,75,167,476,522
Paramnames: N,theta,I,Over

Function Correlate_dists(dists : Context[I,RunIndex] ; rankcorrs : num~~
eric array[I,J] ; ~
  I,J : IndexType;~
  RunIndex : optional Index = Run )
Title: Correlate Dists (d, rc, i, j )
Description: Reorders the samples in dists so as to match the desired~~
 rank correlations between distributions as closely as possible.  Ran~~
kCorrs must be positive definite, and the diagonal should contain all~~
 ones.~
~
The result will be distributions having the same margins as the origi~~
nal input, but with rank correlations close to those of the rankcorrs~~
 matrix.
Definition: if not IsSampleEvalMode and Handle(RunIndex)=Handle(Run) T~~
hen~
   dists  {Mid mode}~
Else begin~
  var u := if Handle(RunIndex)=Handle(run) ~
               Then Sample(Gaussian(0,rankcorrs,I,J))~
               Else Random(Gaussian(0,rankcorrs,I,J),Over:RunIndex);~
  var dsort := sortIndex(dists,RunIndex);~
  var urank := UniqueRank(u,RunIndex);~
  dists[RunIndex=dsort[RunIndex=urank]]~
end
Nodelocation: 136,392,1
Nodesize: 100,16
Windstate: 2,301,193,557,477
Paramnames: dists,rankcorrs,I,J,RunIndex

Function Correlate_with( S, ref : Context[RunIndex] ; rc : scalar ; ~
  RunIndex : optional Index = Run )
Title: Correlate With (s, ref, rc )
Description: Reorders the samples of S so that the result is correlate~~
d with the reference sample with a rank correlation close to rankcorr~~
.  ~
~
Example: To generate a logNormal distribution that is highly correlat~~
ed with Ch1, use, e.g.,:  Correlate_With( LogNormal(2,3), Ch1, 0.8 )~
~
Note: This achieves a given unweighted rank correlation.  If you have~~
 a non-default SampleWeighting of points, the weighted rank correlato~~
n may differ.
Definition: if IsSampleEvalMode or Handle(runIndex)<>Handle(Run) Then ~~
begin~
  Index q := 1..2;~
  var u := If Handle(RunIndex)=Handle(Run) ~
               Then binormal( 0, 1, q, rc )~
                Else Random(binormal(0,1,q,rc),Over:RunIndex);~
  var rrank := UniqueRank(ref,RunIndex);~
  var u1sort := sortIndex(u[q=1],RunIndex);~
  var u2rank := UniqueRank(u[q=2],RunIndex);~
  var ssort := sortIndex(S,RunIndex);~
  S[RunIndex=ssort[RunIndex=u2rank[RunIndex=~
                                                            u1sort[Ru~~
nIndex=rrank]]]]~
end ~
else {mid mode}~
  S
Nodelocation: 128,312,1
Nodesize: 96,16
Windstate: 2,205,170,545,485
Paramnames: S,ref,rc,RunIndex

Function Uniformspherical(I : IndexType ; R : optional Numeric[I,Run] ~~
;~
Over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Uniform Spherical (i, r )
Description: Generates points uniformly on a sphere (or circle or hype~~
rsphere).~
Each sample generated is indexed by I -- so if I has 3 elements, the ~~
points will lie on a sphere.~
~
The mid value is a bit strange here since there isn't really a median~~
 that lies on the sphere.  Obviously the center of the sphere is the ~~
middle value, but that isn't in the allowable range.  So, an arbitrar~~
y point on the sphere is used.
Definition: if IsNotSpecified(R) then R:=1;~
var u := Normal(0,1,over:I,~
                          singleSampleMethod:singleSampleMethod); ~
var d := sqrt( sum(u^2,I) );~
ifall d=0 and @I then R/sqrt(size(I)) else r*u/d
Nodelocation: 328,168,1
Nodesize: 86,16
Windstate: 2,151,227,476,424
Paramnames: I,R,Over

Function Multiuniform(corr : Numeric[I,J,Run] ; I,J : IndexType ; lb,u~~
b : optional Numeric[I,J,Run] ;~
Over : ... optional atomic ;~
singleSampleMethod : hidden optional scalar )
Title: MultiUniform ( c, i, j, lb, ub )
Description: The multi-variate uniform distribution.~
Generates vector samples (indexed by I) such that each component has ~~
a uniform marginal distribution, and such that each component have th~~
e pair-wise correlations given by corr.  Indexes I and J must have th~~
e same number of elements, corr needs to be symmetric and must obey a~~
 certain semidefinite condition (namely that the transformed matrix [~~
 2*sin(30*cov) ] is positive semidefinite.  In most cases, this rough~~
ly the same as corr being, or not being, positive semidefinite).  Lb ~~
and ub can be used to specify upper and lower bounds, either for all ~~
components, or individually if these bounds are indexed by I.  If lb ~~
& ub are omitted, each component will have marginal Uniform(0,1).~
~
The correlation specified in corr is true sample correlation - not ra~~
nk correlation.   ~
~
The transformation here is based on:~
* Falk, M. (1999), "A simple approach to the generation of uniformly ~~
distributed random variables with prescribed correlations," Comm. in ~~
Stats - Simulation and Computation 28: 785-791.
Definition: if IsNotSpecified(lb) then lb:=0;~
if IsNotSpecified(ub) then ub := 1;~
var R := if I=J then 1 else 2*sin(30*corr);~
var g := Gaussian(0,R,I,J,~
                              singleSampleMethod:singleSampleMethod);~~
~
Cumnormal( g ) * (ub-lb) + lb
Nodelocation: 132,168,1
Nodesize: 100,16
Windstate: 2,67,106,608,611
Paramnames: corr,I,J,lb,ub,Over

Module Depricated_multi_var
Title: Depricated multi-variate stuff
Description: Functions found in this module are here for legacy reason~~
s.  They existed in older versions of the Multivariate library, but h~~
ave been become obsolete for whatever reason.
Author: Lonnie
Date: Mon, Apr 30, 2007 3:49 PM
Defaultsize: 48,24
Nodelocation: 80,944,1
Nodesize: 56,32

Function Samplecovariance(X ; I : Index ; J : optional Index ; R : Ind~~
ex)
Title: Sample Covariance
Description: This function is obsolete.  In Analytica 4.0, the builtin~~
 function Variance can be used to compute a covariance matrix.  The e~~
quivalent of this function would be:  Variance( X, R, CoVarDim:I, CoV~~
arDim2:J ).~
~
Returns a covariance matrix based on the sampled data, X, indexed by ~~
I and R.  (I is the dimensionality of X, R corresponds to the samples~~
).  The result will be indexed by I and J -- supply J to be the same ~~
length as I.~
~
Note that the mean is simply Average(X,R), and doen't warrant a separ~~
ate function.
Definition: var I2 := if IsNotSpecified(J) ~
              Then (Index K/((identifier of I)&"2") := I do VarTerm(K~~
)) ~
              Else VarTerm(J);~
var Z:=X-Average(X,R);~
var Zt := Z[@I=@I2];~
Sum(Z*Zt,R)/(size(R)-1)
Nodelocation: 80,48,1
Nodesize: 48,24
Windstate: 2,222,299,476,297
Paramnames: X,I,J,R

Function Samplecorrelation(X : array[I,R] ; I,J,R : IndexType)
Title: sample correlation
Description: This function is obsolete.  A covariance matrix can be co~~
mputed in Analytica 4.0+ using the built-in function Correlation.  Th~~
e equivalent of this function is Correlation(X,X[@I=@J],R).~
~
Returns a correlation matrix based on data in X, where each data poin~~
t is a vector indexed by I, and the entries in the correlation matrix~~
 are the pair-wise correlations of the columns of data.  A second ind~~
ex, J, of size identical to I, is required in order to index the 2-di~~
mensional result.
Definition: var z:=x-average(x,R);~
var zt := slice(z,I,cumulate(1,J));~
sum(z*zt,R) / sqrt(sum(z^2,R) * sum(zt^2,R))~

Nodelocation: 208,48,1
Nodesize: 48,24
Windstate: 2,70,24,523,377
Paramnames: X,I,J,R

Close Depricated_multi_var

Text Multvar_te1
Description: Parametric Multivariate  Distributions
Nodelocation: 160,40,-1
Nodesize: 136,12

Text Multvar_te2
Description: Creating an array of mutually correlated distributions:
Nodelocation: 232,368,-1
Nodesize: 200,16

Text Multvar_te3
Description: Creating a single univariate  distribution correlated wit~~
h another existing dist:
Nodelocation: 296,280,-1
Nodesize: 268,12

Function Normal_correl(m, s, r, y: Numeric ;~
over : optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Normal_correl(m, s, r, y)
Description: Generates a normal distribution with mean m, standard dev~~
iation s, and correlation r with normally distributed value y.  In a ~~
deterministic context, it will return m.~
~
If y is not normally distributed, the result will also not be normal,~~
 and the correlation will be approximate. It generalizes appropriatel~~
y if any of the parameters are arrays:The result array will have the ~~
union of  the indexes of the parameters.
Definition: IF r<-1 OR r>1 THEN Error('Correlation parameter r in func~~
tion Normal_correl(m, s, r, y) is outside the expected range [-1, 1].~~
');~
IFOnly IsSampleEvalMode ~
THEN m + s * (Sqrt(1-r^2) ~
                        * Normal(Sameindexes( 0, m ), Sameindexes( 1,~~
 s ),~
                                       singleSampleMethod:singleSampl~~
eMethod )  ~
                     +  r * (y - Mean(y))/Sdeviation(y))~
ELSE m
Nodelocation: 352,312,1
Nodesize: 108,16
Windstate: 2,102,90,503,416
Paramnames: m,s,r,y,over

Module Multivariate_interna
Title: Multivariate Internal Functions
Author: Lonnie
Date: Tue, May 01, 2007 9:29 PM
Defaultsize: 48,24
Nodelocation: 200,944,1
Nodesize: 52,32
Diagstate: 1,605,145,550,300,17

Function Sameindexes(x, y)
Title: SameIndexes(x,y)
Description: Returns an array with the same indexes as y, and value x ~~
in each cell.
Definition: IF y=y THEN x ELSE x
Nodelocation: 120,64,1
Nodesize: 80,20
Paramnames: x,y

Function Uniquerank(X : Array[I]; I : Index)
Title: UniqueRank
Description: Returns the Rank of X along I, but such that the rank ass~~
igned is unique for every element.  Thus, when there are ties, instea~~
d of getting the same rank, as would happen with the Rank(X,I) functi~~
on, the ranks will be assigned arbitrarily.   Consider:~
   [ 3, 1, 3, 2, 3, 2, 1 ]~
Ranks become:~
   [5,1,6,3,7,4,2 ]
Definition: index Pos := @I;~
var s := SortIndex(X[@I=@Pos],Pos);~
var result := 1;~
for n:=Pos do ( result[@I=s[@Pos=n]] := n );~
result
Nodelocation: 272,64,1
Nodesize: 52,20
Windstate: 2,477,347,537,379
Paramnames: X,I

Close Multivariate_interna

Function Multinormal(m, s: Numeric; cm: ArrayType[i, j,Run]; i , j: In~~
dexType ;~
Over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Multinormal(m,s,c,i,j)
Description: A multi-variate normal (or Gaussian) distribution with me~~
an  m, standard deviation s, and correlation matrix cm.  m and s may ~~
be scalar or indexed by i. cm must be symmetric, positive-definite, a~~
nd indexed by i & j, which must be the same length.~
~
Multinormal uses a correlation matrix.  Compare with Gaussian, which ~~
also defines a multi-variate normal but which uses a covariance matri~~
x.
Definition: Gaussian(m,cm*s*s[@i=@j],i,j,over,singleSampleMethod)
Nodelocation: 472,72,1
Nodesize: 84,16
Windstate: 2,391,248,512,343
Paramnames: m,s,cm,i,j,Over

Text Multvar_te4
Description: Reshaped distributions:
Nodelocation: 136,448,-1
Nodesize: 100,16

Function Dist_reshape(x : Numeric[R] ; newdist : all Numeric[R] ; ~
R : optional Index = Run )
Title: Dist_reshape(x, newdist)
Description: Reshapes the probability distribution of uncertain quanti~~
ty x so that it has the same marginal probability distribution (i.e, ~~
same set of sample values) as newdist, but retains the same ranks as ~~
x.  Thus:~
     Rank(Sample(x), Run) ~
           = Rank(Sample(Reshape_dist(x, y)), Run)~
In a Mid context, it simply returns the mid value of newdist, with an~~
y indexes of x.~
~
The result retains any rank correlations that x may have with other p~~
redecessor variables.  So, the rank-order correlation between a third~~
 variable z and x will be the same as the rank-order correlation betw~~
een z and a reshaped version of x, i.e.~
     RankCorrel(x, z) = RankCorrel(Reshape_Dist(x, y), z)~
~
The operation may optionally be applied along an index other than Run~~
.
Definition: IFOnly IsSampleEvalMode or Handle(R)<>Handle(Run) THEN BEG~~
IN~
     VAR dsort := SortIndex(newdist, Run);~
     VAR xranks := Rank(x, Run);~
     newdist[Run = dsort[Run=xranks]]~
     END~
ELSE newdist * (x=x)
Nodelocation: 152,472,1
Nodesize: 116,16
Windstate: 2,102,90,646,469
Paramnames: x,newdist,R

Text Multvar_te5
Description: Arrays with serial correlation
Nodelocation: 208,532,-1
Nodesize: 168,12

Function Normal_serial_correl(m, s, r: Numeric; i: IndexType ;~
over : ... optional atomic;~
singleSampleMethod : optional hidden scalar )
Title: Normal_serial_correl(m,s,r,i)
Description: Generates an array over index i of normal distributions w~~
ith mean m, standard deviation s, and correlation r between successiv~~
e values  over index i.  You can give each distribution a different m~~
ean and/or standard deviation if m and/or s are arrays indexed by i. ~~
If r is indexed by i, r[i=k] specifies the correlation between result~~
[i=k] and result[i=k-1].  (Then the first correlation, slice(r, i, 1)~~
  is ignored.)
Definition: Var x := Normal(0, 1,singleSampleMethod:singleSampleMethod~~
);~
(FOR j := i DO  ~
       x := Normal_correl( 0, 1, r[i = j],x,~
                                       singleSampleMethod:singleSampl~~
eMethod ) ) ~
* s + m
Nodelocation: 160,560,1
Nodesize: 120,16
Windstate: 2,353,325,540,383
Paramnames: m,s,r,i,over

Function Normal_additive_gro(x, m, s, r: Numeric; i: IndexType ;~
over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Normal_additive_gro(x,m,s,r,i)
Description: Adds a normally distributed percent growth g with mean m ~~
and standard deviation s to x for each value of index i. The growth g~~
 for each i has serial correlation r with g for i-1.
Definition: x *( 1 + Cumulate(Normal_serial_correl(m, s, r, i,~
                            singleSampleMethod:singleSampleMethod), i~~
))
Nodelocation: 159,600,1
Nodesize: 119,16
Windstate: 2,102,90,519,306
Paramnames: x,m,s,r,i,over

Function Normal_compound_gro(x, m, s, r: Numeric; t: IndexType ;~
over : ... optional atomic;~
singleSampleMethod : optional hidden scalar )
Title: Normal_compound_gro(x,m,s,r,t)
Description: An array of values over time index t, starting from with ~~
value x, and with compound growth applied for each time interval, wit~~
h normal uncertainty with mean m and standard deviation s  The growth~~
 g for each i has correlation r with g for i-1.
Definition: x * Cumproduct(IF t = Slice(t, 1) THEN 1 ELSE Normal_seria~~
l_correl(m, s, r, t, singleSampleMethod:singleSampleMethod ) + 1, t)~~

Nodelocation: 159,640,1
Nodesize: 119,16
Windstate: 2,102,90,529,366
Paramnames: x,m,s,r,t,over

Function Dist_serial_correl(x; r; i: IndexType ;~
over : ... optional atomic;~
singleSampleMethod : optional hidden scalar )
Title: Dist_serial_correl(x,r,i)
Description: Generates an array y over index i where each y[i] has a m~~
arginal  distribution identical to x, and serial rank correlation of ~~
r with y[i-1].  If x  is indexed by i,  each y[i] has the same margin~~
al distribution as x[i], but with samples reordered to have the speci~~
fied rank correlation r between successive values.  If r is indexed b~~
y i, r[i=k] specifies the rank correlation between y[i=k] and y[i=k-1~~
].  Then the first correlation, r[i=1], is ignored.~
~
In Mid context, it returns Mid(x).~
~
Note: The result retains no probabilistic dependence on x.
Definition: Dist_reshape(Normal_serial_correl( 0, 1, r, i, singleSampl~~
eMethod:singleSampleMethod ), x)
Nodelocation: 408,560,1
Nodesize: 120,16
Windstate: 2,302,78,477,447
Paramnames: x,r,i,over

Function Dist_additive_growth(x, g, r: Numeric; i: IndexType;~
over : ... optional atomic;~
singleSampleMethod : optional hidden scalar )
Title: Dist_additive_growth(x,g,r,i)
Description: Generates an array of values over index i, with the first~~
 equal to x, and successive values adding an uncertain growth with pr~~
obability distribution g, and serial correlation r between growth[i =~~
 k] and growth[i=k-1].  x, g, and r each may be indexed by i if you w~~
ant them to vary over i.
Definition: x + Cumulate(Dist_serial_correl( g, r, i, singleSampleMeth~~
od : singleSampleMethod), i)
Nodelocation: 407,600,1
Nodesize: 119,16
Windstate: 2,102,90,506,300
Paramnames: x,g,r,i,over

Function Dist_compound_growth(x, g, r; i: IndexType ;~
over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Dist_compound_growth(x,g,r,i)
Description: Starts with x and applies a compound growth g for each va~~
lue of index i. The growth g for each i has correlation r with g for ~~
i-1.
Definition: x * Cumproduct(~
           IF i = Slice(i, 1) THEN 1 ~
           ELSE (Dist_serial_correl( g, r, i,   ~
                                    singleSampleMethod:singleSampleMe~~
thod ) + 1)~
, i)
Nodelocation: 407,640,1
Nodesize: 119,16
Windstate: 2,102,90,489,307
Paramnames: x,g,r,i,over

Text Multvar_te6
Description: Distributions on Linear Regression coefficients
Nodelocation: 296,688,-1
Nodesize: 256,12

Function Regressionnoise( Y : Numeric[I,Run] ; B : Numeric[I,K,Run] ; ~~
I,K : Index; C : optional Numeric[K,Run] )
Title: RegressionNoise(Y,B,I,K,C)
Description: When you have data, Y[I] and B[I,K], generated from an un~~
derlying model with unknown coefficients C[k] and S of the form:~
~
Y = Sum( C*B, I) + Normal(0,S)~
~
This function computes an estimate for S.  ~
~
When using in conjunction with RegressionDist, it is most efficient t~~
o provide the optional parameter C to both routines, where C is the e~~
xpected value of the regression coefficients, obtained from calling R~~
egression(Y,B,I,K).  Doing so avoids an unnecessary call to the built~~
in Regression function.
Definition: if IsNotSpecified(C) Then C := Regression(Y,B,I,K);~
Var resid := Y - Sum(C*B,K);~
sqrt( Sum(resid^2,I) / (size(I)-size(K)) );~

Nodelocation: 384,736,1
Nodesize: 104,20
Windstate: 2,332,211,498,542
Paramnames: Y,B,I,K,C

Function Regressionfitprob( Y : Numeric[I,Run] ; B : Numeric[I,K,Run] ~~
; I,K : Index; C : optional Numeric[K,Run] ; ~
S :  optional Numeric[I,Run] )
Title: RegressionFitProb(Y,B,I,K,C)
Description: Once you've obtained regression coefficients C (indexed b~~
y K) by calling the Regression function, this function returns the pr~~
obability that a fit this poor would occur by chance, given the assum~~
ption that the data was generated by a process of the form:~
~
 Y = Sum( C*B,K) + Normal(0,S)~
~
If this result is very close to zero, it probably indicates that the ~~
assumption of linearity is bad.  If it is very close to one, then it ~~
validates the assumption of linearity.~
~
This is not a distribution function - it does not return a sample whe~~
n evaluated in Sample mode.  However, it does complement the multivar~~
iate RegressionDist function also included in this library.~
~
To use, first call the Regression function, then you must either know~~
 the measurement knows a priori, or obtain it using the RegressionNoi~~
se function.~
~
Var E_C := Regression(Y,B,I,K);~
Var S := RegressionNoise(Y,B,I,K,C);~
Var PrThisPoor := RegressionFitProb(Y,B,I,K,E_C,S)
Definition: if IsNotSpecified(C) then C:=Regression(Y,B,I,K);~
if IsNotSpecified(S) then S:=RegressionNoise(Y,B,I,K);~
var resid := Y - sum(C*B,K);~
var n := size(I);~
var chi2 := sum( resid^2 / Mean(S)^2, I);~
GammaI( n/2 - 1, chi2/2 )
Nodelocation: 152,800,1
Nodesize: 112,20
Windstate: 2,287,69,586,548
Paramnames: Y,B,I,K,C,S

Close Multivariate_distrib