Updating GitHub Pages

CBOR/docs/PeterO.Cbor.CBORException.html

@@ -8,8 +8,8 @@ <h2>PeterO.Cbor.CBORException</h2>
 
 <pre>public class CBORException :
     System.Exception,
-    System.Runtime.InteropServices._Exception,
-    System.Runtime.Serialization.ISerializable
+    System.Runtime.Serialization.ISerializable,
+    System.Runtime.InteropServices._Exception
 </pre>
 
 <p>Exception thrown for errors involving CBOR data.</p>

CBOR/docs/PeterO.Cbor.CBORType.html

@@ -8,9 +8,9 @@ <h2>PeterO.Cbor.CBORType</h2>
 
 <pre>public sealed struct CBORType :
     System.Enum,
+    System.IConvertible,
     System.IFormattable,
-    System.IComparable,
-    System.IConvertible
+    System.IComparable
 </pre>
 
 <p>Represents a type that a CBOR object can have.</p>

CBOR/docs/PeterO.Rounding.html

@@ -8,9 +8,9 @@ <h2>PeterO.Rounding</h2>
 
 <pre>public sealed struct Rounding :
     System.Enum,
+    System.IConvertible,
     System.IFormattable,
-    System.IComparable,
-    System.IConvertible
+    System.IComparable
 </pre>
 
 <p><b>Deprecated.</b> Use ERounding from PeterO.Numbers/com.upokecenter.numbers.</p>

CBOR/docs/PeterO.TrapException.html

@@ -8,8 +8,8 @@ <h2>PeterO.TrapException</h2>
 
 <pre>public class TrapException :
     System.ArithmeticException,
-    System.Runtime.InteropServices._Exception,
-    System.Runtime.Serialization.ISerializable
+    System.Runtime.Serialization.ISerializable,
+    System.Runtime.InteropServices._Exception
 </pre>
 
 <p><b>Deprecated.</b> Use ETrapException from PeterO.Numbers/com.upokecenter.numbers.</p>

MailLib/docs/PeterO.Mail.MessageDataException.html

@@ -8,8 +8,8 @@ <h2>PeterO.Mail.MessageDataException</h2>
 
 <pre>public class MessageDataException :
     System.Exception,
-    System.Runtime.InteropServices._Exception,
-    System.Runtime.Serialization.ISerializable
+    System.Runtime.Serialization.ISerializable,
+    System.Runtime.InteropServices._Exception
 </pre>
 
 <p>Exception thrown when a message has invalid syntax.</p>

MailLib/docs/PeterO.Text.Normalization.html

@@ -8,9 +8,9 @@ <h2>PeterO.Text.Normalization</h2>
 
 <pre>public sealed struct Normalization :
     System.Enum,
+    System.IConvertible,
     System.IFormattable,
-    System.IComparable,
-    System.IConvertible
+    System.IComparable
 </pre>
 
 <p>Represents a Unicode normalization form.</p>

Numbers/docs/PeterO.Numbers.ERounding.html

@@ -8,9 +8,9 @@ <h2>PeterO.Numbers.ERounding</h2>
 
 <pre>public sealed struct ERounding :
     System.Enum,
+    System.IConvertible,
     System.IFormattable,
-    System.IComparable,
-    System.IConvertible
+    System.IComparable
 </pre>
 
 <p>Specifies the mode to use when &quot;shortening&quot; numbers that otherwise can&#39;t fit a given number of digits, so that the shortened number has about the same value. This &quot;shortening&quot; is known as rounding. (The &quot;E&quot; stands for &quot;extended&quot;, and has this prefix to group it with the other classes common to this library, particularly EDecimal, EFloat, and ERational.).</p>

Numbers/docs/PeterO.Numbers.ETrapException.html

@@ -8,8 +8,8 @@ <h2>PeterO.Numbers.ETrapException</h2>
 
 <pre>public sealed class ETrapException :
     System.ArithmeticException,
-    System.Runtime.InteropServices._Exception,
-    System.Runtime.Serialization.ISerializable
+    System.Runtime.Serialization.ISerializable,
+    System.Runtime.InteropServices._Exception
 </pre>
 
 <p>Exception thrown for arithmetic trap errors. (The &quot;E&quot; stands for &quot;extended&quot;, and has this prefix to group it with the other classes common to this library, particularly EDecimal, EFloat, and ERational.).</p>

randomfunc.html

@@ -1279,7 +1279,7 @@ <h3>Correlated Random Numbers</h3>
 <li>calculate <code>[x, y*sqrt(1 - rho * rho) + rho * x]</code>, where <code>rho</code> is a <em>correlation coefficient</em> in the interval [-1, 1] (if <code>rho</code> is 0, the variables are uncorrelated).</li>
 </ul>
 
-<p>Another way to generate correlated random numbers is explained in the section &quot;<a href="#Gaussian_and_Other_Copulas"><strong>Gaussian and Other Copulas</strong></a>&quot;.</p>
+<p>Other ways to generate correlated random numbers are explained in the section &quot;<a href="#Gaussian_and_Other_Copulas"><strong>Gaussian and Other Copulas</strong></a>&quot;.</p>
 
 <p><a id=Specific_Non_Uniform_Distributions></a></p>
 
@@ -1899,7 +1899,7 @@ <h3>Gaussian and Other Copulas</h3>
 <li><code>erf(v)</code> is the <a href="https://en.wikipedia.org/wiki/Error_function"><strong>error function</strong></a> of the variable <code>v</code> (see the appendix).</li>
 </ul>
 
-<p>&mdash;</p>
+<p>&nbsp;</p>
 
 <pre>METHOD GaussianCopula(covar)
    mvn=MultivariateNormal(nothing, covar)
@@ -1921,8 +1921,8 @@ <h3>Gaussian and Other Copulas</h3>
 <p>Other kinds of copulas describe different kinds of correlation between random numbers.  Examples of other copulas are&mdash;</p>
 
 <ul>
-<li>the <strong>Fr&eacute;&ndash;Hoeffding upper bound copula</strong> <em>[x, x, ..., x]</em> (e.g., <code>[x, x]</code>), where <code>x = RNDU01()</code>,</li>
-<li>the <strong>Fr&eacute;&ndash;Hoeffding lower bound copula</strong> <code>[x, 1.0 - x]</code> where <code>x = RNDU01()</code>,</li>
+<li>the <strong>Fr&eacute;chet&ndash;Hoeffding upper bound copula</strong> <em>[x, x, ..., x]</em> (e.g., <code>[x, x]</code>), where <code>x = RNDU01()</code>,</li>
+<li>the <strong>Fr&eacute;chet&ndash;Hoeffding lower bound copula</strong> <code>[x, 1.0 - x]</code> where <code>x = RNDU01()</code>,</li>
 <li>the <strong>product copula</strong>, where each number is a separately generated <code>RNDU01()</code> (indicating no correlation between the numbers), and</li>
 <li>the <strong>Archimedean copulas</strong>, described by M. Hofert and M. M&auml;chler (2011)<sup><a href="#Note13"><strong>(13)</strong></a></sup>.</li>
 </ul>
@@ -1969,7 +1969,7 @@ <h3>Other Non-Uniform Distributions</h3>
 <li><strong>Chi distribution</strong>: <code>sqrt(GammaDist(df * 0.5, 2))</code>, where <code>df</code> is the number of degrees of freedom.</li>
 <li><strong>Cosine distribution</strong>: <code>min + (max - min) * atan2(x, sqrt(1 - x * x)) / pi</code>, where <code>x = RNDNUMRANGE(-1, 1)</code> and <code>min</code> is the minimum value and <code>max</code> is the maximum value (Saucier 2000, p. 17; inverse sine replaced with <code>atan2</code> equivalent).</li>
 <li><strong>Double logarithmic distribution</strong>: <code>min + (max - min) * (0.5 + (RNDINT(1) * 2 - 1) * 0.5 * RNDU01OneExc() * RNDU01OneExc())</code>, where <code>min</code> is the minimum value and <code>max</code> is the maximum value (see also Saucier 2000, p. 15, which shows the wrong X axes).</li>
-<li><strong>Frech&eacute;t distribution</strong>: <code>b*pow(-ln(RNDU01ZeroExc()),-1.0/a) + loc</code>, where <code>a</code> is the shape, <code>b</code> is the scale, and <code>loc</code> is the location of the distribution&#39;s curve peak (mode). This expresses a distribution of maximum values.</li>
+<li><strong>Fr&eacute;chet distribution</strong>: <code>b*pow(-ln(RNDU01ZeroExc()),-1.0/a) + loc</code>, where <code>a</code> is the shape, <code>b</code> is the scale, and <code>loc</code> is the location of the distribution&#39;s curve peak (mode). This expresses a distribution of maximum values.</li>
 <li><strong>Generalized extreme value (Fisher&ndash;Tippett) distribution</strong>: <code>a - (pow(-ln(RNDU01ZeroOneExc()), -c) - 1) * b / c</code> if <code>c != 0</code>, or <code>a - ln(-ln(RNDU01ZeroOneExc())) * b</code> otherwise, where <code>b</code> is the scale, <code>a</code> is the location of the distribution&#39;s curve peak (mode), and <code>c</code> is a shape parameter. This expresses a distribution of maximum values.</li>
 <li><strong>Generalized Tukey lambda distribution</strong>: <code>(s1 * (pow(x, lamda1)-1.0)/lamda1 - s2 * (pow(1.0-x, lamda2)-1.0)/lamda2) + loc</code>, where <code>x</code> is <code>RNDU01()</code>, <code>lamda1</code> and <code>lamda2</code> are shape parameters, <code>s1</code> and <code>s2</code> are scale parameters, and <code>loc</code> is a location parameter.</li>
 <li><strong>Half-normal distribution</strong>: <code>abs(Normal(0, sqrt(pi * 0.5) / invscale)))</code>, where <code>invscale</code> is a parameter of the half-normal distribution.</li>
@@ -1992,7 +1992,7 @@ <h3>Other Non-Uniform Distributions</h3>
 <li><strong>Zeta distribution</strong>: Generate <code>n = floor(pow(RNDU01ZeroOneExc(), -1.0 / r))</code>, and if <code>d / pow(2, r) &lt; (d - 1) * RNDU01OneExc() * n / (pow(2, r) - 1.0)</code>, where <code>d = pow((1.0 / n) + 1, r)</code>, repeat this process. The parameter <code>r</code> is greater than 0. Based on method described in Devroye 1986. A zeta distribution <a href="#Censoring_and_Truncation"><strong>truncated</strong></a> by rejecting random values greater than some positive integer is called a <em>Zipf distribution</em> or <em>Estoup distribution</em>. (Note that Devroye uses &quot;Zipf distribution&quot; to refer to the untruncated zeta distribution.)</li>
 </ul>
 
-<p>The Python sample code also contains implementations of the <strong>power normal distribution</strong>, the <strong>power lognormal distribution</strong>, the <strong>negative multinomial distribution</strong>,  the <strong>multivariate <em>t</em>-distribution</strong>, and the <strong>multivariate Poisson distribution</strong>.</p>
+<p>The Python sample code also contains implementations of the <strong>power normal distribution</strong>, the <strong>power lognormal distribution</strong>, the <strong>negative multinomial distribution</strong>,  the <strong>multivariate <em>t</em>-distribution</strong>, the <strong>multivariate <em>t</em>-copula</strong>, and the <strong>multivariate Poisson distribution</strong>.</p>
 
 <p><a id=Geometric_Sampling></a></p>
 

randomfunc.md

@@ -1033,7 +1033,7 @@ According to (Saucier 2000), sec. 3.8, to generate two correlated (dependent) ra
 - generate two independent and identically distributed random variables `x` and `y` (for example, two `Normal(0, 1)` variables or two `RNDU01()` variables), and
 - calculate `[x, y*sqrt(1 - rho * rho) + rho * x]`, where `rho` is a _correlation coefficient_ in the interval \[-1, 1\] (if `rho` is 0, the variables are uncorrelated).
 
-Another way to generate correlated random numbers is explained in the section "[**Gaussian and Other Copulas**](#Gaussian_and_Other_Copulas)".
+Other ways to generate correlated random numbers are explained in the section "[**Gaussian and Other Copulas**](#Gaussian_and_Other_Copulas)".
 
 <a id=Specific_Non_Uniform_Distributions></a>
 ## Specific Non-Uniform Distributions
@@ -1575,7 +1575,7 @@ One example is a _Gaussian copula_; this copula is sampled by sampling from a [*
 - The parameter `covar` is the covariance matrix for the multinormal distribution.
 - `erf(v)` is the [**error function**](https://en.wikipedia.org/wiki/Error_function) of the variable `v` (see the appendix).
 
-&mdash;
+&nbsp;
 
     METHOD GaussianCopula(covar)
        mvn=MultivariateNormal(nothing, covar)
@@ -1593,8 +1593,8 @@ Each of the resulting uniform numbers will be in the interval [0, 1], and each o
 
 Other kinds of copulas describe different kinds of correlation between random numbers.  Examples of other copulas are&mdash;
 
-- the **Fr&eacute;&ndash;Hoeffding upper bound copula** _\[x, x, ..., x\]_ (e.g., `[x, x]`), where `x = RNDU01()`,
-- the **Fr&eacute;&ndash;Hoeffding lower bound copula** `[x, 1.0 - x]` where `x = RNDU01()`,
+- the **Fr&eacute;chet&ndash;Hoeffding upper bound copula** _\[x, x, ..., x\]_ (e.g., `[x, x]`), where `x = RNDU01()`,
+- the **Fr&eacute;chet&ndash;Hoeffding lower bound copula** `[x, 1.0 - x]` where `x = RNDU01()`,
 - the **product copula**, where each number is a separately generated `RNDU01()` (indicating no correlation between the numbers), and
 - the **Archimedean copulas**, described by M. Hofert and M. M&auml;chler (2011)<sup>[**(13)**](#Note13)</sup>.
 
@@ -1636,7 +1636,7 @@ Miscellaneous:
 - **Chi distribution**: `sqrt(GammaDist(df * 0.5, 2))`, where `df` is the number of degrees of freedom.
 - **Cosine distribution**: `min + (max - min) * atan2(x, sqrt(1 - x * x)) / pi`, where `x = RNDNUMRANGE(-1, 1)` and `min` is the minimum value and `max` is the maximum value (Saucier 2000, p. 17; inverse sine replaced with `atan2` equivalent).
 - **Double logarithmic distribution**: `min + (max - min) * (0.5 + (RNDINT(1) * 2 - 1) * 0.5 * RNDU01OneExc() * RNDU01OneExc())`, where `min` is the minimum value and `max` is the maximum value (see also Saucier 2000, p. 15, which shows the wrong X axes).
-- **Frech&eacute;t distribution**: `b*pow(-ln(RNDU01ZeroExc()),-1.0/a) + loc`, where `a` is the shape, `b` is the scale, and `loc` is the location of the distribution's curve peak (mode). This expresses a distribution of maximum values.
+- **Fr&eacute;chet distribution**: `b*pow(-ln(RNDU01ZeroExc()),-1.0/a) + loc`, where `a` is the shape, `b` is the scale, and `loc` is the location of the distribution's curve peak (mode). This expresses a distribution of maximum values.
 - **Generalized extreme value (Fisher&ndash;Tippett) distribution**: `a - (pow(-ln(RNDU01ZeroOneExc()), -c) - 1) * b / c` if `c != 0`, or `a - ln(-ln(RNDU01ZeroOneExc())) * b` otherwise, where `b` is the scale, `a` is the location of the distribution's curve peak (mode), and `c` is a shape parameter. This expresses a distribution of maximum values.
 - **Generalized Tukey lambda distribution**: `(s1 * (pow(x, lamda1)-1.0)/lamda1 - s2 * (pow(1.0-x, lamda2)-1.0)/lamda2) + loc`, where `x` is `RNDU01()`, `lamda1` and `lamda2` are shape parameters, `s1` and `s2` are scale parameters, and `loc` is a location parameter.
 - **Half-normal distribution**: `abs(Normal(0, sqrt(pi * 0.5) / invscale)))`, where `invscale` is a parameter of the half-normal distribution.
@@ -1658,7 +1658,7 @@ Miscellaneous:
 - **Tukey lambda distribution**: `(pow(x, lamda)-pow(1.0-x,lamda))/lamda`, where `x` is `RNDU01()` and `lamda` is a shape parameter (if 0, the result is a logistic distribution).
 - **Zeta distribution**: Generate `n = floor(pow(RNDU01ZeroOneExc(), -1.0 / r))`, and if `d / pow(2, r) < (d - 1) * RNDU01OneExc() * n / (pow(2, r) - 1.0)`, where `d = pow((1.0 / n) + 1, r)`, repeat this process. The parameter `r` is greater than 0. Based on method described in Devroye 1986. A zeta distribution [**truncated**](#Censoring_and_Truncation) by rejecting random values greater than some positive integer is called a _Zipf distribution_ or _Estoup distribution_. (Note that Devroye uses "Zipf distribution" to refer to the untruncated zeta distribution.)
 
-The Python sample code also contains implementations of the **power normal distribution**, the **power lognormal distribution**, the **negative multinomial distribution**,  the **multivariate _t_-distribution**, and the **multivariate Poisson distribution**.
+The Python sample code also contains implementations of the **power normal distribution**, the **power lognormal distribution**, the **negative multinomial distribution**,  the **multivariate _t_-distribution**, the **multivariate _t_-copula**, and the **multivariate Poisson distribution**.
 
 <a id=Geometric_Sampling></a>
 ## Geometric Sampling

randomgen.py

@@ -627,6 +627,61 @@ def multivariate_t(self, mu, cov, df):
      cd=self.gamma(df*0.5,2.0/df)
      return [(0 if mu==None else mu[i])+mn[i]/math.sqrt(cd) for i in range(len(mn))]
 
+  def _pochhammer(self, a, b):
+    return math.gamma(a+b)/math.gamma(a)
+
+  def _beta(self, a, b):
+    return math.gamma(a)*math.gamma(b)/math.gamma(a+b)
+
+  def _betainc(self, x, a, b):
+    # Incomplete beta function.  NOTE: The SciPy method
+    # scipy.stats.betainc(a, b, x) is the same as _betainc(x, a, b).
+    if x>0.5 and x < 1.0: return 1.0 - betainc(1.0 - x, b, a)
+    if x==0 and a>0: return 0.0
+    if b<50 and math.floor(b)==b:
+        if b<0: return 0
+        return (x**a)*sum([ \
+            self._pochhammer(a,i)*pow(1-x,i)*1.0/math.gamma(i+1) \
+            for i in range(int(b))])
+    if a>0 and a<50 and math.floor(a)==a:
+        return 1.0-((1.0-x)**b)*sum([ \
+            self._pochhammer(b,i)*(x**i)*1.0/math.gamma(i+1) \
+            for i in range(int(a))])
+    ret=pow(10,-100)
+    d=0
+    c=ret
+    i=0
+    k=0
+    while i < 100:
+        # Get next convergent of continued fraction
+        if i==0: num=1.0
+        else:
+          if (i&1)==1: num=-(a+k)*(a+b+k)*x*1.0/((a+i-1)*(a+i))
+          else: num=(b-k)*k*x*1.0/((a+i-1)*(a+i))
+        c=1+num/c # 1 is the convergent's denominator
+        d=1+num*d # ditto
+        if d==0: d=pow(10,-100)
+        if c==0: c=pow(10,-100)
+        d=1.0/d
+        delta=d*c
+        ret*=delta
+        if abs(delta-1.0)<pow(10,-14): break
+        i=i+1
+        if (i&1)==0: k=k+1
+    return ret*(x**a)*((1-x)**b)/(a*self._beta(a,b))
+
+  def _student_t_cdf(self,nu, x):
+    if x<=0:
+       return self._betainc(nu/(x*x+nu),nu*0.5,0.5)*0.5
+    else:
+       return (self._betainc((x*x)/(x*x+nu),0.5,nu*0.5)+1)*0.5
+
+  def t_copula(self, cov, df):
+    """ Multivariate t-copula. 'cov' is the covariance matrix
+       and 'df' is the degrees of freedom.  """
+    mt=self.multivariate_t(None, cov, df)
+    return [self._student_t_cdf(df, c) for c in mt]
+
   def randomwalk_u01(self,n):
      """ Random walk of uniform 0-1 random numbers. """
      ret=[0 for i in range(n+1)]