HP17bII+ Programming t-distribution

[daisyclover18]

Hi, does anyone know how to program the student's t-distribution on this calculator's solver?

[Don Shepherd]

(07-23-2015 06:25 AM)daisyclover18 Wrote:  Hi, does anyone know how to program the student's t-distribution on this calculator's solver?

Daisy, according to the formula given at Wikipedia, this should work, but doesn't, at least on by 17bii:

(fact(ip(((v+1)/2)-1))

/

(sqrt(pi x v) x fact(ip((v/2)-1))))

x

(((1+sq(t)/v)) ^ (-((v+1)/2))))

that's a lot of parentheses, but if they didn't match it wouldn't "compile" and it does, and the error occurs when you CALC, enter a value for v and then solve for t. It says "solution not found" pretty quickly. Perhaps someone can look at this and see why it does not work.

[daisyclover18]

Thanks Don. Hope someone can look into this.

[Don Shepherd]

I tried reducing the number of parentheses, keeping in mind the order of operations, and the equation is now a bit simpler, but it still does not work, you get the same error message.

I'm going to play around with it some more today, there has got to be something I am missing here.

fact(ip((v+1)/2)-1)

/

(sqrt(pi x v) x fact(ip(v/2)-1))

x

((1+t x t/v)^-((v+1)/2))

[SlideRule]

(07-23-2015 06:25 AM)daisyclover18 Wrote:  Hi, does anyone know how to program the student's t-distribution on this calculator's solver?

If I understand the question correctly, you are seeking the equation for

[attachment=2363]

my interpretation follows

T = (XB1-XB2) / SQRT( (SQ(S1)/N1) + (SQ(S2)/N2) )

BEST!

SlideRule

[Thomas Klemm]

Code:

T=IF(MOD(N:2)=0:
  FACT(N-1)÷2^(N-1)÷SQ(FACT(N÷2-1)):
  2^(N-1)×SQ(FACT((N-1)÷2))÷FACT(N-1)÷PI)
÷SQRT(N)
÷(1+SQ(X)÷N)^((N+1)÷2)

Examples:
PDF[StudentTDistribution[7],1.6]
N=7
X=1.6
⇒ T=0.1106654

PDF[StudentTDistribution[4],0.3]
N=4
X=0.3
⇒ T=0.3547096

Kind regards
Thomas

[Dieter]

(07-23-2015 06:25 AM)daisyclover18 Wrote:  Hi, does anyone know how to program the student's t-distribution on this calculator's solver?

What exactly are you looking for?
The PDF? The CDF? The quantile function? One-sided or two-sided?

There can be no helpful solution unless these questions have been answered.

Dieter

[Dieter]

(07-23-2015 09:52 PM)Don Shepherd Wrote:  Daisy, according to the formula given at Wikipedia, this should work,

No, it shouldn't. ;-)

(07-23-2015 09:52 PM)Don Shepherd Wrote:  but doesn't,

Of course it doesn't.

The problem is in the Gamma function. Yes, for integer (!) arguments Γ(x) can be replaced by the factorial of x–1. But Student's t-distribution always deals with half-integers, either in the nominator or in the denominator. You always have to evaluate something like Γ(2.5)/Γ(2) or Γ(7)/Γ(6.5). That's why you cannot simply remove the fractional part by using IP.

As far as I can see the 17BII does not offer a Gamma function (unlie many other HPs). But there is a way: Γ(x+0.5) can be calculated from Γ(x), and there even is a simple way to evaluate Γ(x+0.5)/Γ(x). Thomas' solution uses such an approach (that's the part with "2^(N-1)").

On calculators without Gamma function (e.g. the '41) I like the following method for evaluating the quotient g = Γ(n+0.5)/Γ(n):

Code:

Let [n] = IP(n)
Let  q  = (2[n]-1)! / 2^(2[n]-1) / ([n]-1)!^2 * sqrt(pi)

If n = [n]
   THEN g = q
   ELSE g = [n]/q

Please note that in our case (Student's distribution) n equals half the number of degrees of freedom, i.e. n = ν/2.

For very large n (that may cause overflow) another approach can be useful: a few terms of a simple power series will do.

Dieter

[Thomas Klemm]

(07-23-2015 09:52 PM)Don Shepherd Wrote:  according to the formula given at Wikipedia

\(\frac{\Gamma \left(\frac{\nu+1}{2} \right)} {\sqrt{\nu\pi}\,\Gamma \left(\frac{\nu}{2} \right)} \left(1+\frac{x^2}{\nu} \right)^{-\frac{\nu+1}{2}}\)

We can use the following identities:

• \(\Gamma(n)=(n-1)!\) for \(n\epsilon \mathbb{N}\)
• \(\Gamma(x+1)=x \Gamma(x)\)
• \(\Gamma\left(\tfrac{1}{2}\right)=\sqrt{\pi}\)
  

For \(\nu\) even:

\(\begin{align}
\Gamma(\frac{\nu+1}{2})&=\frac{\nu-1}{2}\cdot\frac{\nu-3}{2}\cdots\frac{1}{2}\cdot\Gamma\left(\tfrac{1}{2}\right) \\
&=\frac{\nu-1}{2}\cdot\frac{\nu-2}{\nu-2}\cdot\frac{\nu-3}{2}\cdot\frac{\nu-4}{\nu-4}\cdots\frac{2}{2}\cdot\frac{1}{2}\cdot\Gamma\left(\tfrac{1}{2}\right) \\
&=\frac{\nu-1}{2}\cdot\frac{\nu-2}{2\cdot(\frac{\nu}{2}-1)}\cdot\frac{\nu-3}{2}\cdot\frac{\nu-4}{2\cdot(\frac{\nu}{2}-2)}\cdots\frac{2}{2\cdot1}\cdot\frac{1}{2}\cdot\Gamma\left(\tfrac{1}{2}\right) \\
&=\frac{(\nu-1)!}{2^{\nu-1}\,(\frac{\nu}{2}-1)!}\,\sqrt{\pi} \\
\end{align}\)

\(\Gamma(\frac{\nu}{2})=(\frac{\nu}{2}-1)!\)

\(\begin{align}
\frac{\Gamma(\frac{\nu+1}{2})} {\sqrt{\pi}\,\Gamma(\frac{\nu}{2})}&=\frac{(\nu-1)!}{\sqrt{\pi}\,2^{n-1}\,(\frac{\nu}{2}-1)!\,(\frac{\nu}{2}-1)!}\,\sqrt{\pi} \\
&=\frac{(\nu-1)!}{2^{n-1}\,((\frac{\nu}{2}-1)!)^2} \\
\end{align}\)

For \(\nu\) odd:

\(\Gamma(\tfrac{\nu+1}{2})=(\tfrac{\nu+1}{2}-1)!=\frac{\nu-1}{2}!\)

\(\begin{align}
\Gamma(\frac{\nu}{2})&=\frac{\nu-2}{2}\cdot\frac{\nu-4}{2}\cdots\frac{1}{2}\cdot\Gamma\left(\tfrac{1}{2}\right) \\
&=\frac{\nu-1}{\nu-1}\cdot\frac{\nu-2}{2}\cdot\frac{\nu-3}{\nu-3}\cdot\frac{\nu-4}{2}\cdots\frac{2}{2}\cdot\frac{1}{2}\cdot\Gamma\left(\tfrac{1}{2}\right) \\
&=\frac{\nu-1}{2\cdot\frac{\nu-1}{2}}\cdot\frac{\nu-2}{2}\cdot\frac{\nu-3}{2\cdot(\frac{\nu-1}{2}-1)}\cdot\frac{\nu-4}{2}\cdots\frac{2}{2\cdot1}\cdot\frac{1}{2}\cdot\Gamma\left(\tfrac{1}{2}\right)​ \\
&=\frac{(\nu-1)!}{2^{\nu-1}\,\frac{\nu-1}{2}!}\,\sqrt{\pi} \\
\end{align}\)

\(\begin{align}
\frac{\Gamma(\frac{\nu+1}{2})} {\sqrt{\pi}\,\Gamma(\frac{\nu}{2})}&=\frac{\frac{\nu-1}{2}!\,2^{\nu-1}\,\frac{\nu-1}{2}!}{\sqrt{\pi}\,(\nu-1)!\,\sqrt{\pi}} \\
&=\frac{2^{\nu-1}\,(\frac{\nu-1}{2}!)^2}{(\nu-1)!\,\pi} \\
\end{align}\)

Kind regards
Thomas

[Thomas Klemm]

(07-25-2015 08:06 AM)Dieter Wrote:  For very large n (that may cause overflow) another approach can be useful: a few terms of a simple power series will do.

Or then just use the normal distribution.

[Don Shepherd]

Thanks Thomas. I think I am getting outside of my area of comfort!

I tried just taking that original Wikipedia formula and translating it directly to 17b solverese. When I said "it didn't work," I meant that it resulted in "solution not found" when I tried to run it (but I gather from yours and Dieter's info that it wouldn't have worked even if it "worked"). I played around with it, like changing the negative subscript to 1/ the positive subscript, and a couple of other things, but I never could escape the dreaded "solution not found." Something in my implementation of that equation causes the 17b to not be able to evaluate it at all. Oh well.

Thanks again.

[Thomas Klemm]

(07-25-2015 11:29 AM)Don Shepherd Wrote:  I think I am getting outside of my area of comfort!

You might still remember Valentin's HP-15C Mini-challenge that was solved by Gerson using Euler's reflection formula:

\(\Gamma(z)\Gamma(1-z) = \frac{\pi}{\sin{(\pi z)}}\)

Now if you set \(z=\frac{1}{2}\) you get:

\(\Gamma(\tfrac{1}{2})\Gamma(\tfrac{1}{2}) = \frac{\pi}{\sin{(\tfrac{\pi}{2})}} = \frac{\pi}{1} = \pi \)

From this we can conclude that in fact: \(\Gamma(\tfrac{1}{2})=\sqrt{\pi}\)

Often I have to write down concrete examples to understand formulas.
Thus for \(\nu=6\) we can use \(\Gamma(x+1)=x \Gamma(x)\) to get:

\(\begin{align}
\Gamma(\frac{\nu+1}{2})&=\Gamma(\frac{7}{2})=\Gamma(\frac{5}{2}+1) \\
&=\frac{5}{2}\cdot\Gamma(\frac{5}{2}) \\
&=\frac{5}{2}\cdot\frac{3}{2}\cdot\Gamma(\frac{3}{2}) \\
&=\frac{5}{2}\cdot\frac{3}{2}\cdot\frac{1}{2}\cdot\Gamma(\frac{1}{2}) \\
&=\frac{5}{2}\cdot\frac{3}{2}\cdot\frac{1}{2}\cdot\sqrt{\pi} \\
\end{align}\)

We don't end up with a factorial but with the double factorial 5!! in the numerator.
Thus we just insert the missing values which happen to be even:

\(\begin{align}
\Gamma(\frac{7}{2})&=\frac{5}{2}\cdot\frac{3}{2}\cdot\frac{1}{2}\cdot\sqrt{\pi} \\
&=\frac{5}{2}\cdot\frac{4}{4}\cdot\frac{3}{2}\cdot\frac{2}{2}\cdot\frac{1}{2}​\cdot\sqrt{\pi} \\
&=\frac{5}{2}\cdot\frac{4}{2\cdot2}\cdot\frac{3}{2}\cdot\frac{2}{2\cdot 1}\cdot\frac{1}{2}\cdot\sqrt{\pi} \\\end{align}\)

Now we end up with 5! in the nominator. In the denominator we end up with the product of 2⁵ and 2!:

\(\Gamma(\frac{7}{2})=\frac{5!}{2^5\,2!}\,\sqrt{\pi}\)

Compare this to the general formula if you set \(\nu=6\):

\(\Gamma(\frac{6+1}{2})=\frac{(6-1)!}{2^{6-1}\,(\frac{6}{2}-1)!}\,\sqrt{\pi}\)

Quote:Something in my implementation of that equation causes the 17b to not be able to evaluate it at all.

I didn't enter your formula but as it's written right now the T=(…) is missing.
I assume that in this case the expression is just set to 0.
And this might explain that no solution was found.

HTH
Thomas

[Don Shepherd]

Yes, if you omit the = the solver simply equates what you give it to 0 and solves for the unknown variable. So it should solve for t in my example, but it gives that error message.

I'm out of my league in this thread. I'm going to start a new thread with a 17b formula for calculating a standard mod-10 check digit, like they do on credit card numbers and other ID numbers. I wrote a solver equation for that years ago and recently improved on it.

[Thomas Klemm]

(07-25-2015 05:21 PM)Don Shepherd Wrote:  So it should solve for t in my example, but it gives that error message.

Oops, I didn't notice that you are using the variable t.
What you are trying to calculate is the value of a function f(v, t)=(…).
But this function has no roots t where f(v, t)=0 as it's strictly positive.
Thus you'd have to use a different variable e.g. F=(…).

Cheers
Thomas

[Don Shepherd]

(07-25-2015 05:54 PM)Thomas Klemm Wrote:  

(07-25-2015 05:21 PM)Don Shepherd Wrote:  So it should solve for t in my example, but it gives that error message.

Oops, I didn't notice that you are using the variable t.
What you are trying to calculate is the value of a function f(v, t)=(…).
But this function has no roots t where f(v, t)=0 as it's strictly positive.
Thus you'd have to use a different variable e.g. F=(…).

Cheers
Thomas

Yeah, I see what you mean. I just didn't know enough about this student t distribution thing. I saw the formula in Wikipedia and it seemed like the only variables were v (degrees of freedom) and t, and I read somewhere that you specify the degrees of freedom, therefore I thought you need to solve for t, but now I see that you specify a value for t also (so there are really two inputs) and then the value you are looking for is the whole expression with t and v plugged in. So X = (the whole expression) and solve for x.

But I've already removed my equation from my 17b after I saw your work and Dieter's work.

[rprosperi]

(07-25-2015 05:21 PM)Don Shepherd Wrote:  I'm out of my league in this thread...

Me too, sheesh!!

I've followed it too, as the Pioneer solver always interests me, but when Dieter and Thomas kicked-in, I feel like I'm standing in jet-wash behind an F-15 taking off.

You guys must be post-grad Math Professors, right? Please don't tell me "no, it's just a hobby" as I'll be forever crushed...

I don't follow it all, but I am always interested and impressed with your contributions and level of commitment (just to type some of those replies is impressive). And despite my ignorance in these areas, I feel like I am always learning a little, even if only my osmosis.

[Thomas Klemm]

(07-25-2015 09:24 PM)rprosperi Wrote:  You guys must be post-grad Math Professors, right?

It's true that I have a mathematical background. But I'm using math mostly for recreational reasons. There's still a lot of fascinating stuff that I don't know and I'm eager to learn.

Quote:I feel like I am always learning a little, even if only my osmosis.

That's good to know as it makes me a little sad that we've lost both of you.
And probably daisyclover18 as well:

Quote:Last Visit: 24.07.2015 13:45

Kind regards
Thomas

[rprosperi]

(07-28-2015 06:52 PM)Thomas Klemm Wrote:  It's true that I have a mathematical background. But I'm using math mostly for recreational reasons. There's still a lot of fascinating stuff that I don't know and I'm eager to learn.

Thomas - As an engineer, I'd always used and enjoyed mathematics as a fundamental tool for just about all aspects of any engineering analysis, but couldn't see any point in spending endless hours doing math just for math's sake. After spending time with Joe Horn at PPC meetings, he introduced me to basic aspects of number theory; mostly as reasons to write programs, but eventually for it's own sake. So while I certainly have no background in much of the math you guys do here, I've at least opened my mind enough to be curious about it, and every now and then, I learn something. So bay all means, please continue, even if it feels like you may be losing the crowd. You never know when the silent audience is learning...