{ From user Lonnie, Model Delaunay_from_python at 1-Aug-2018 5:32:51 PM, encoding="UTF-8" }
SoftwareVersion 5.2.0


{ System Variables with non-default values: }
SampleSize := 1000
TypeChecking := 1
Checking := 1
SaveOptions := 2
SaveValues := 0

NodeInfo FormNode: 1,0,0,1,0,0,0,,0,0,,0,0



Model Delaunay_from_python
Description: This example demonstrates how an Analytica model can call methods in Python (Python is a programming language) using COM. This example is from the artilce~
~
• <a href="http://wiki.analytica.com/index.php?title=COM_Integration/Calling_Python_code">Calling Python code</a> ~
~
that appears on the Analytica wiki. In this example, a Delaunay tessellation, aka Delaunay triangulation, is computed for a random set of 2-D points. I Delaunay tessellation is a set of adjacent and non-overlapping triangles, using the points as vertices, such that any circle that contains the points of a triangle does not contain any other point in its interior. 
Author: Lonnie Chrisman, Ph.D.~
Lumina Decision Systems
Date: Wed, Aug 1, 2018 9:25 AM
SaveAuthor: Lonnie
SaveDate: Wed, Aug 1, 2018 5:32 PM
DiagState: 2,1,0,773,545,17
WindState: 2,412,462,720,350
FontStyle: Arial,15
FileInfo: 0,Model Delaunay_from_python,2,2,0,0,C:\Users\Lonnie\Documents\Analytica\Analytica_Python\Delaunay from python.ana

Variable py
Title: py
Definition: COMCreateObject("Lumina.DelaunayCOM")
NodeLocation: 160,80,1
NodeSize: 48,24

Button Pause
Title: Pause
NodeLocation: 160,176,1
NodeSize: 48,24
OnClick: py->Pause()

Index Pt
Title: Pt
Description: The point number. Indexes the data points.
Definition: 1..10
NodeLocation: 312,80,1
NodeSize: 48,24
WindState: 2,389,524,720,350

Index Dim
Title: Dim
Description: Points in this example are 2-D. Indexes the point dimensions. I.e., Dim=1 is the point's first dimension, Dim=2 is the 2nd dimension.
Definition: [1,2]
NodeLocation: 424,80,1
NodeSize: 48,24
WindState: 2,408,498,720,350
Att_PrevIndexValue: [1,2]

Variable pts
Title: pts
Description: A random set of 2-D points that will be used.
Definition: Random( Binormal( Array(Dim,[1,2]), Array(Dim,[2,1]), Dim, 0.5, over:Pt ) )
NodeLocation: 312,168,1
NodeSize: 48,24
WindState: 2,501,586,720,350
ValueState: 2,6,315,634,576,1,MIDM
GraphSetup: Att_ContLineStyle Graph_Primary_Valdim:4
Att_CoordinateIndex: Dim

{!40404|FreePassObjectCount 86sQq6qJq4RNUM0tPBWqCVgcLOsQ2PZu9SyrmK$cu3RVdU8zkXgw6NLRRBxazLn8X2yvd9g9RsJsqnfU8WzaSmryodN$Tj2X7DGDtYyZrNSWehdciIrISdXeQSKuXuTRccTNEwc9Ugz1vfUTjPsLZt_wjZM_les9OOCunLoWwMZomwvTC$9LdsxpZPGiJhr662FpP5kDVfv6D1sQBhBa5t5452aH1aoJiMSJC1R0e39ZqMK3qXpGm31KTueLt7Sww030gLAGo5QUfPK0U9V5JWjbaPxo0GckqvkUR$hQYt1DNB6gLtaYy7KGB1qSoDoewz23fQEr7lCnlnkiRtRkq31GmiP8p3UyD9SZygLr3Mo0wxjHwi0Qn20$8ZJ0Yo7bZUJ70eNlNbu2JskZDb2em_473$kDzJk0PwqohbHf8OQXnZRG$RuVEMYfnP_TxUm9BTLI6yf4jUes$9rdHuOuBEYSRH9uL0yKTelRICvH1aTTPQ8_kEk7QegmYb_UwRf$FFF8pR8hAWZwYcRIubtKqjhzBEYIpMe1$8KvZYWFFhu75AqfXo8pUvbjcXByIpNKMky3RFqR_6PluipwtnXCoJl_UQ7sXnCgxt7GPRMD$cBYAE3wj5eEdh4SkvvqgT7cMYVUPvaIp$X1QfrvwndIpbpopmI$lIW5d_JVdfYP6fTjknmK3rQgIsJcs_3_udC2KNST3qer7iHi4KVYYTDqi05CFtgWl3hJmBThoqjaD7TajoSH9QmR0VzKYklmbGCaju$hYSl7rY2Yt7HLJAtrFQdmUNJe2pR2VvBQZWM5$7ThphL61vmS1Y0Pfk50OGtEb_NoDf6a3Y3b7gFqQ1fIxcI_hP9vfRE2siYQIB50yusrqqsux$4AGOWgq0CPcr5MdwDXtCZxKk9a0UyRxS_X5fFrT5kO3lR9tdN9xkYND4xqjeZVSQPPQRUXchnu09JTfs3HWn1JcwEawIg3TuKnFjDjEnKuT2fHvZDvbI0lVG2qeSH7_slfaVSPNNNOQSWbgnu09KVht5KZq5Og$Kg1OnBb0SvNsMtPxV3fFrT6mR6pWF_lWH4tiYOG81xtpmjiiijmpty29HPZjv5IWl$GXq8Sn7UsFf4VyQuNuPxU2cCoQ2hM1jQ8scM8wjXMD4xqkeaVTQPPQSVYdiow2BKVhu5KZp5NfzHd_Lk7XzOrJoHoKsP_Z8mN0gK$iP7scN9xlZPF6$uniebYXWWXadhlry3CLWhs3HWm1IbuCXtEc$OpEh8c5b6eAkIuU6kN1iO4oXG0nZMB$riZSLGB742111368DIOVdmw4GSft6MdvBVq9WtFe2TvLpHnHoLuS1cDrT7nS8rZI2oaMA_pfVNG94$yvustuvy$39FMUeoz8LZn1GXp6Ql5SpCb0RtJmGlGnKtS1dEsU7nS8rYH0nZMA_qgYQJD83$ywvuvx_27DJRaju3FSft6MdvCWrBYwJj7Z0TxQwQxT0a9lM_dHyeL3nXI3rdQF4wnfYSNJGEDDEGIMQVahow3EOan0FVm2LfzJf1OnAb0SvNsMtPzX6iIwZCtXDvdL4qbNA$qgXPIC73$ywvuuwy$38EMUdny8LZo1IYq7Qk4Qn9YyNqHlEkFnJtQ$bBoR4kP5pWG0nZMA$qfWOHA50yv}
Library 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 functions in this library.  A Mid value from a distribution function will therefore be indexed by I, whlie a Sample from a distribution function is indexed by both I and Run.  These distribution functions can also 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 distributions.  Correlate_with, for example, allows you to generate a univarite 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, where 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: Wed, May 10, 2017 3:31 PM
DefaultSize: 48,24
NodeLocation: 640,64,1
NodeSize: 64,24
NodeInfo: 1,1,1,1,1,1,0,0,0,0,,,,0
DiagState: 2,15,10,743,821,17
WindState: 2,401,199,776,387
FontStyle: Arial, 15
Att_PreLoadScript: {!40404|FreePassObjectCount 86sQq6qJq4RNUM0tPBWqCVgcLOsQ2PZu9SyrmK$cu3RVdU8zkXgw6NLRRBxazLn8X2yvd9g9RsJsqnfU8WzaSmryodN$Tj2X7DGDtYyZrNSWehdciIrISdXeQSKuXuTRccTNEwc9Ugz1vfUTjPsLZt_wjZM_les9OOCunLoWwMZomwvTC$9LdsxpZPGiJhr662FpP5kDVfv6D1sQBhBa5t5452aH1aoJiMSJC1R0e39ZqMK3qXpGm31KTueLt7Sww030gLAGo5QUfPK0U9V5JWjbaPxo0GckqvkUR$hQYt1DNB6gLtaYy7KGB1qSoDoewz23fQEr7lCnlnkiRtRkq31GmiP8p3UyD9SZygLr3Mo0wxjHwi0Qn20$8ZJ0Yo7bZUJ70eNlNbu2JskZDb2em_473$kDzJk0PwqohbHf8OQXnZRG$RuVEMYfnP_TxUm9BTLI6yf4jUes$9rdHuOuBEYSRH9uL0yKTelRICvH1aTTPQ8_kEk7QegmYb_UwRf$FFF8pR8hAWZwYcRIubtKqjhzBEYIpMe1$8KvZYWFFhu75AqfXo8pUvbjcXByIpNKMky3RFqR_6PluipwtnXCoJl_UQ7sXnCgxt7GPRMD$cBYAE3wj5eEdh4SkvvqgT7cMYVUPvaIp$X1QfrvwndIpbpopmI$lIW5d_JVdfYP6fTjknmK3rQgIsJcs_3_udC2KNST3qer7iHi4KVYYTDqi05CFtgWl3hJmBThoqjaD7TajoSH9QmR0VzKYklmbGCaju$hYSl7rY2Yt7HLJAtrFQdmUNJe2pR2VvBQZWM5$7ThphL61vmS1Y0Pfk50OGtEb_NoDf6a3Y3b7gFqQ1fIxcI_hP9vfRE2siYQIB50yusrqqsux$4AGOWgq0CPcr5MdwDXtCZxKk9a0UyRxS_X5fFrT5kO3lR9tdN9xkYND4xqjeZVSQPPQRUXchnu09JTfs3HWn1JcwEawIg3TuKnFjDjEnKuT2fHvZDvbI0lVG2qeSH7_slfaVSPNNNOQSWbgnu09KVht5KZq5Og$Kg1OnBb0SvNsMtPxV3fFrT6mR6pWF_lWH4tiYOG81xtpmjiiijmpty29HPZjv5IWl$GXq8Sn7UsFf4VyQuNuPxU2cCoQ2hM1jQ8scM8wjXMD4xqkeaVTQPPQSVYdiow2BKVhu5KZp5NfzHd_Lk7XzOrJoHoKsP_Z8mN0gK$iP7scN9xlZPF6$uniebYXWWXadhlry3CLWhs3HWm1IbuCXtEc$OpEh8c5b6eAkIuU6kN1iO4oXG0nZMB$riZSLGB742111368DIOVdmw4GSft6MdvBVq9WtFe2TvLpHnHoLuS1cDrT7nS8rZI2oaMA_pfVNG94$yvustuvy$39FMUeoz8LZn1GXp6Ql5SpCb0RtJmGlGnKtS1dEsU7nS8rYH0nZMA_qgYQJD83$ywvuvx_27DJRaju3FSft6MdvCWrBYwJj7Z0TxQwQxT0a9lM_dHyeL3nXI3rdQF4wnfYSNJGEDDEGIMQVahow3EOan0FVm2LfzJf1OnAb0SvNsMtPzX6iIwZCtXDvdL4qbNA$qgXPIC73$ywvuuwy$38EMUdny8LZo1IYq7Qk4Qn9YyNqHlEkFnJtQ$bBoR4kP5pWG0nZMA$qfWOHA50yv}
{!40400|Att_clearTypeFonts: -1}

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) distribution, 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 matrix.  ~
~
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 ignoring sampling error.  That may be insignificant of N is large.  Otherwise, 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, after 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:singleSampleMethod))~
             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

Function Invertedwishart(psi:Number[I,J,Run] ; n :positive ; I,J : Index ;~
singleSampleMethod : optional hidden scalar)
Title: Inverted Wishart(Psi,n,I,J)
Description: The inverted Wishart distribution represents a distribution of covariance matrices.  Each sample from the InvertedWishart is a covariance matrix.  It is conjugate to a Wishart distribution, and therefore is sometimes used as a Bayesian prior distribution for covariance.  The parameter, Psi, must be a positive definite matrix.~
~
Suppose you represent the prior distribution of covariance using an inverted Wishart distribution: InvertedWishart(Psi,m).  You observe some data, X[I,R], where R:=1..N indexes each datapoint and I is the vector 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 about 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, singleSampleMethod:singleSampleMethod), I, J )
NodeLocation: 352,224,1
NodeSize: 108,16
WindState: 2,102,90,553,530

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]).  You 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 itself uncertain.  This function returns the coefficients C as a distribution (i.e., in Sample mode, it returns a sampling of coefficients indexed by Run and K), reflecting the uncertainty in the estimation of these 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 follows:~
~
{ 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 parameters 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(RegressionNoise(...)) 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 parameter.
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 Linv := Decompose(Sum( B*B[K=K2]/S^2,I),K,K2);~
var L := Invert(LInv,K,K2);~
var Z := Normal(0,1,over:K2,~
				singleSampleMethod:singleSampleMethod);~
Sum( L*Z,K2 ) + C
NodeLocation: 152,736,1
NodeSize: 112,20
WindState: 2,947,-10,561,882

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 whose elements sum to 1.  It is commonly used to represent second order probability 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 confidence 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: 296,120,1
NodeSize: 64,16
WindState: 2,26,18,624,485

Function Binormal(MeanVec :numeric[I,Run]; Sdeviations : positive[I,Run]; 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 the indicated individual standard deviations (>0) and the indicated correlation coefficient.  The index, I, must have exactly 2 elements, Sdeviations 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: 88,16
WindState: 2,2,24,525,540

Function Multinomial(N:NonNegative ; theta:NonNegative  ; I : IndexType;~
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 distribution 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 the corresponding outcome (die number) occurred during that sample point.  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: 120,120,1
NodeSize: 96,16
WindState: 2,75,167,476,522

Function Correlate_dists(dists : Context[I,RunIndex] ; rankcorrs : numeric 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.  RankCorrs must be positive definite, and the diagonal should contain all ones.~
~
The result will be distributions having the same margins as the original input, but with rank correlations close to those of the rankcorrs matrix.
Definition: if not IsSampleEvalMode and Handle(RunIndex)=Handle(Run) Then~
   dists  {Mid mode}~
Else begin~
  var corr := if @I=@J then 1 else 2 * Sin(30 * rankCorrs);~
  var u := if Handle(RunIndex)=Handle(run) ~
               Then Sample(Gaussian(0,corr,I,J))~
               Else Random(Gaussian(0,corr ,I,J),Over:RunIndex);~
  var dsort := sortIndex(dists,RunIndex,position:true);~
  var urank := Rank(u,RunIndex,type:null);~
  dists[@RunIndex=dsort[@RunIndex=urank]]~
end
NodeLocation: 136,400,1
NodeSize: 104,16
WindState: 2,693,59,721,615

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 correlated with the reference sample with a rank correlation close to rankcorr.  ~
~
Example: To generate a logNormal distribution that is highly correlated 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 correlaton may differ.
Definition: if IsSampleEvalMode or Handle(runIndex)<>Handle(Run) Then begin~
  Index q := 1..2;~
  var corr := 2 * sin( 30 * rc );~
  var u := If Handle(RunIndex)=Handle(Run) ~
               Then binormal( 0, 1, q, corr )~
                Else Random(binormal(0,1,q, corr),Over:RunIndex);~
  var rrank := Rank(ref,RunIndex,type:Null);~
  var u1sort := sortIndex(u[q=1],RunIndex,position:true);~
  var u2rank := Rank(u[q=2],RunIndex,type:Null);~
  var ssort := sortIndex(S,RunIndex,position:true);~
  S[@RunIndex=ssort[@RunIndex=u2rank[@RunIndex=~
                                                            u1sort[@RunIndex=rrank]]]]~
end ~
else {mid mode}~
  S
NodeLocation: 144,312,1
NodeSize: 112,16
WindState: 2,816,36,715,607

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 hypersphere).~
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 arbitrary 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: 352,168,1
NodeSize: 96,16
WindState: 2,151,227,476,424

Function Multiuniform(corr : Numeric[I,J,Run] ; I,J : IndexType ; lb,ub : 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 the pair-wise correlations given by corr.  Indexes I and J must have the 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 roughly 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 rank 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: 136,168,1
NodeSize: 112,16
WindState: 2,67,106,608,611

Module Depricated_multi_var
Title: Depricated multi-variate stuff
Description: Functions found in this module are here for legacy reasons.  They existed in older versions of the Multivariate library, but have 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 : Index)
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 equivalent of this function would be:  Variance( X, R, CoVarDim:I, CoVarDim2: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 separate 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

Function Samplecorrelation(X : array[I,R] ; I,J,R : IndexType)
Title: sample correlation
Description: This function is obsolete.  A covariance matrix can be computed in Analytica 4.0+ using the built-in function Correlation.  The equivalent of this function is Correlation(X,X[@I=@J],R).~
~
Returns a correlation matrix based on data in X, where each data point 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 index, J, of size identical to I, is required in order to index the 2-dimensional 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

Close Depricated_multi_var

Text Multvar_te1
Description: Parametric Multivariate  Distributions
NodeLocation: 220,42,-1
NodeSize: 196,14

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 with another existing dist:
NodeLocation: 314,278,-1
NodeSize: 286,10

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 deviation 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 appropriately 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 function 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:singleSampleMethod )  ~
                     +  r * (y - Mean(y))/Sdeviation(y))~
ELSE m
NodeLocation: 368,312,1
NodeSize: 108,16
WindState: 2,102,90,503,416

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: 56,32
DiagState: 1,605,145,518,268,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

Function Uniquerank(X : Array[I]; I : Index)
Title: UniqueRank
Description: Returns the Rank of X along I, but such that the rank assigned is unique for every element.  Thus, when there are ties, instead of getting the same rank, as would happen with the Rank(X,I) function, the ranks will be assigned arbitrarily.   Consider:~
   [ 3, 1, 3, 2, 3, 2, 1 ]~
Ranks become:~
   [5,1,6,3,7,4,2 ]
Definition: Rank(X,I,type:Null)
NodeLocation: 272,64,1
NodeSize: 56,20
WindState: 2,477,347,537,379

Close Multivariate_interna

Function Multinormal(m, s: Numeric; cm: ArrayType[i, j,Run]; i , j: IndexType ;~
Over : ... optional atomic ;~
singleSampleMethod : optional hidden scalar )
Title: Multinormal(m,s,c,i,j)
Description: A multi-variate normal (or Gaussian) distribution with mean  m, standard deviation s, and correlation matrix cm.  m and s may be scalar or indexed by i. cm must be symmetric, positive-definite, and 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 matrix.
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

Text Multvar_te4
Description: Reshaped distributions:
NodeLocation: 136,448,-1
NodeSize: 100,16

Function Dist_reshape(x : Numeric[R] ; newdist : Numeric[R] ; ~
R : optional Index = Run )
Title: Dist_reshape(x, newdist)
Description: Reshapes the probability distribution of uncertain quantity 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 any indexes of x.~
~
The result retains any rank correlations that x may have with other predecessor variables.  So, the rank-order correlation between a third variable z and x will be the same as the rank-order correlation between 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 BEGIN~
     VAR dsort := SortIndex(newdist, R, position:true);~
     VAR xranks := Rank(x, R, type:Null);~
     newdist[@R = dsort[@R=xranks]]~
     END~
ELSE newdist * (x=x)
NodeLocation: 152,472,1
NodeSize: 116,16
WindState: 2,102,90,836,582

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 with mean m, standard deviation s, and correlation r between successive values  over index i.  You can give each distribution a different mean 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:singleSampleMethod ) ) ~
* s + m
NodeLocation: 160,560,1
NodeSize: 120,16
WindState: 2,353,325,540,383

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

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, with 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_serial_correl(m, s, r, t, singleSampleMethod:singleSampleMethod ) + 1, t)
NodeLocation: 160,640,1
NodeSize: 128,16
WindState: 2,102,90,529,366

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 marginal  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 marginal distribution as x[i], but with samples reordered to have the specified rank correlation r between successive values.  If r is indexed by 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, singleSampleMethod:singleSampleMethod ), x)
NodeLocation: 440,560,1
NodeSize: 120,16
WindState: 2,302,78,477,447

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 probability 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 want them to vary over i.
Definition: x + Cumulate(Dist_serial_correl( g, r, i, singleSampleMethod : singleSampleMethod), i)
NodeLocation: 439,600,1
NodeSize: 119,16
WindState: 2,102,90,506,300

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 value 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:singleSampleMethod ) + 1)~
, i)
NodeLocation: 439,640,1
NodeSize: 119,16
WindState: 2,102,90,489,307

Text Multvar_te6
Description: Distributions on Linear Regression coefficients
NodeLocation: 296,700,-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 underlying 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 to provide the optional parameter C to both routines, where C is the expected value of the regression coefficients, obtained from calling Regression(Y,B,I,K).  Doing so avoids an unnecessary call to the builtin 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: 112,20
WindState: 2,332,211,498,542

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 by K) by calling the Regression function, this function returns the probability that a fit this poor would occur by chance, given the assumption 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 when evaluated in Sample mode.  However, it does complement the multivariate 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 RegressionNoise 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

Text Multivar_Te24
NodeLocation: 379,256,-1
NodeSize: 380,1
NodeInfo: 1,0,0,0,1,1,0,,0,,0
NodeColor: 0,0,0

Text Multivar_Te38
NodeLocation: 376,341,-1
NodeSize: 380,1
NodeInfo: 1,0,0,0,1,1,0,,0,,0
NodeColor: 0,0,0

Text Multivar_Te39
NodeLocation: 380,424,-1
NodeSize: 380,1
NodeInfo: 1,0,0,0,1,1,0,,0,,0
WindState: 2,102,83,464,224
NodeColor: 0,0,0

Text Multivar_Te40
NodeLocation: 380,505,-1
NodeSize: 380,1
NodeInfo: 1,0,0,0,1,1,0,,0,,0
NodeColor: 0,0,0

Text Multvar_te173623456
Description: Parametric Multivariate  Distributions
NodeLocation: 220,42,-1
NodeSize: 196,14

Text Multvar_te217657837
Description: Creating an array of mutually correlated distributions:
NodeLocation: 232,368,-1
NodeSize: 200,16

Text Multvar_te388723474
Description: Creating a single univariate  distribution correlated with another existing dist:
NodeLocation: 314,278,-1
NodeSize: 286,10

Text Multvar_te4873647383
Description: Reshaped distributions:
NodeLocation: 136,448,-1
NodeSize: 100,16

Text Multvar_te5837466373
Description: Arrays with serial correlation
NodeLocation: 208,532,-1
NodeSize: 168,12

Text Multvar_te6873647383
Description: Distributions on Linear Regression coefficients
NodeLocation: 296,700,-1
NodeSize: 256,12

Text Multvar_Te1736453672
NodeLocation: 379,256,-1
NodeSize: 380,1
NodeInfo: 1,0,0,0,1,1,0,,0,,0
NodeColor: 0,0,0

Text Multvar_Te2837465637
NodeLocation: 376,341,-1
NodeSize: 380,1
NodeInfo: 1,0,0,0,1,1,0,,0,,0
NodeColor: 0,0,0

Text Multvar_Te3338748393
NodeLocation: 380,424,-1
NodeSize: 380,1
NodeInfo: 1,0,0,0,1,1,0,,0,,0
NodeColor: 0,0,0

Text Multvar_Te4928344443
NodeLocation: 380,505,-1
NodeSize: 380,1
NodeInfo: 1,0,0,0,1,1,0,,0,,0
NodeColor: 0,0,0

Text Multvar_Te5273647710
NodeLocation: 380,673,-1
NodeSize: 380,1
NodeInfo: 1,0,0,0,1,1,0,,0,,0
NodeColor: 0,0,0

Alias Gaussian1494358961
Title: Gaussian(m, cv, i, j)
NodeLocation: 104,72,1
NodeSize: 80,16
NodeInfo: 1,1,1,1,1,1,0,0,0,0,,0,1,0
NodeColor: 58978,52427,65535
Original: Gaussian

Close Multivariate_distrib

Variable tessellationVertices
Title: tessellation vertex points
Description: Calls a Python function to compute a Delaunay tessellation of the points. The result value is a list of triangles, where each cell value specifies the point (the value in the Pt index) for the triangle vertex.
Definition: Local tri := py->Tessellation(COMArray(pts,Pt,Dim) )[@.dim2=@Vertex_pt] + 1;~
Triangle := 1..IndexLength(tri.dim1);~
tri[@.dim1=@Triangle]
NodeLocation: 432,168,1
NodeSize: 56,32
WindState: 2,548,541,720,350
ValueState: 2,831,50,421,416,0,MIDM
ReformVal: [Vertex_pt,Triangle]

Index Triangle
Title: Triangle
Description: Indexes the triangles resulting from the tessellation. The value is the triangle number from 1 to number of  triangles.
Definition: ComputedBy(tessellationVertices)
NodeLocation: 312,240,1
NodeSize: 48,24
WindState: 2,376,545,720,350
ValueState: 2,376,379,416,419,,MIDM
{!40300|Att_SlicerPopupSize: 227,210}

Index Vertex_pt
Title: Vertex pt
Description: Indexes the 3 vertices of each triangle.
Definition: 1..3
NodeLocation: 312,296,1
NodeSize: 48,24
WindState: 2,296,515,720,350
{!40300|Att_SlicerPopupSize: 232,80}

Variable closed_tesselation_p
Title: closed tesselation points
Description: The 3 points of each triangle is expanded here to a sequence of 4 points where the first and last points are the same. This facilitates drawing a line plot, since each curve completes the edges of the triangle.
Definition: tessellationVertices[@Vertex_pt=Mod(@Closed_vertex_pt-1,3)+1]
NodeLocation: 432,264,1
NodeSize: 48,32
WindState: 2,410,516,720,350
ValueState: 2,468,473,421,413,,MIDM
ReformVal: [Vertex_pt,Triangle]
Att_ResultSliceState: [Seg,1,Triangle,1,Vertex_pt,1]

Index Closed_vertex_pt
Title: Closed vertex pt
Description: Indexes a path around a triangle, starting and ending at the same vertex. The closed path is useful for drawing each triangle in graph view, since each curve corresponds to 1 triangle.
Definition: 1..4
NodeLocation: 312,352,1
NodeSize: 48,24

Objective Tessellation_plot
Title: Tessellation plot
Description: A line-plot showing the triangles of the tessellation. Note that the edges of adjacent triangles overlap.
Definition: pts[Pt=closed_tesselation_p]
NodeLocation: 608,264,1
NodeSize: 56,24
ValueState: 2,560,351,787,559,1,MIDM
GraphSetup: {!50200|Att_AreaFill Graph_Primary_Valdim:3}
ReformVal: [Closed_vertex_pt,Triangle]
Att_CoordinateIndex: Dim

Close Delaunay_from_python