QuantLabId's profile - activity

I don't think we're going to be able to fix this very easily at all. The problem is Maxima's framework, which there has been much discussion about and which is not going to change.

(%i5) assume(n>0);
(%o5)                               [n > 0]
(%i6) display2d:false;

(%o6) false
(%i7) integrate(sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)),s,0,inf);

Is (n+1)/2 an integer?

no
;
(%o7) gamma(n/2+1/2)*2^(-n/2-1/2)*(t^2/(2*n)+1/2)^(-n/2-1/2)
 /(sqrt(n)*gamma(n/2)*sqrt(pi))
(%i8) integrate(sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)),s,0,inf);

Is (n+1)/2 an integer?

yes;
Is n an integer?

yes;
(%o8) gamma(n/2+1/2)*2^(-n/2-1/2)*(t^2/(2*n)+1/2)^(-n/2-1/2)
 /(sqrt(n)*gamma(n/2)*sqrt(pi))

Don't even ask about why it asks if n is an integer when it knows (n+1)/2 is an integer; Maxima itself acknowledges that its assumption framework is weak, and of course every open source project has limited resources to fix big architectural issues like that. I don't think it would be trivial for them to fix that.

As for Sage, assuming n is an integer will help a little, but only so much.

(%i9) declare(n,integer);

(%o9) done
(%i10) integrate(sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)),s,0,inf);

Is (n+1)/2 an integer?

yes;
(%o10) gamma(n/2+1/2)*2^(-n/2-1/2)*(t^2/(2*n)+1/2)^(-n/2-1/2)
 /(sqrt(n)*gamma(n/2)*sqrt(pi))

(%i11) declare((n+1)/2,integer);

declare: improper argument: (1+n)/2
 -- an error. To debug this try: debugmode(true);

So we wouldn't even be able to try that, even if we wanted to; Maxima only lets "atoms" be declared to be something (i.e., variables), though we can assume positive or negative about them. And even though the answers are the same, many have pointed out that going through a decision tree of this kind and trying to check that can be arbitrarily computationally expensive, if it's even possible.

I'm sorry, but this is something at the core of how Maxima does business. On the plus side, the answer works nicely enough that one gets from Maxima, and you could use it in Sage.

The problem is I want the integration for general n(positive integer), not for some specific n=3 or 5!

Thank you for the answer. But the full solution should be f(t,n). we could consider f(t,n) from f(t,1), f(t,2), f(t,3)... but it's not the exact solution!!!

Yes, it works! But what about when n is an even integer? It also has solution when n is an even number, which is of course p.d.f of t-distribution!

Here's a code for student t pdf.

var("s n t")
assume(n, 'integer')
assume(n>0)
tdist=integral(sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)),s,0,infinity)
tdist

And it says, TypeError: Computation failed since Maxima requested additional constraints (try the command 'assume(n+1)/2>0)' before integral or limit evaluation, for example): Is (n+1)/2 an integer?

The funny thing is when I add n=5 or n=3, there's no error!

var("s n t")
assume(n, 'integer')
n=3
tdist=integral(sqrt(s/n)*1/sqrt(2*pi)*exp(-s*t^2/(2*n))*exp(-s/2)*s^(n/2-1)/(gamma(n/2)*2^(n/2)),s,0,infinity)
tdist

Could any one help me about this?

EDIT : Actually, I tried this on MATHLAB to fail.

   My last hope would be Mathematica? wolframalpha.com says it's timed out!

Firstly, I defined normal distribution density function

sage: def n_d(mean, variance, x): return 1/(sqrt(2pivariance))exp(-1/2(x-mean)^2/variance)

and I tried to double integrate the joint density function n_d*n_d

sage: integrate(integrate(n_d(mean_x,var_x,x)n_d(beta_1x+beta_0,var_e, y),x,-oo,oo),y,-oo,oo)

first n_d is normal distribution of x second n_d is conditional normal distribution on x, which is normal distribution of y

and it came out as not 1

the result was

sage: 1/2sqrt(2beta_1^2var_x + 2var_e)sqrt(2)/sqrt(beta_1^2var_x +var_e)

What could possibly be wrong?

Note that the type of output depends on the base ring of the matrix, for example a floating point (double) matrix

sage: matrix(RDF,[[1/2,-1/4],[-1/2,1/3]]).eigenvalues()
[0.779908245295, 0.0534250880383]

yields a floating-point result. A matrix over the rationals

sage: matrix(QQ,[[1/2,-1/4],[-1/2,1/3]]).eigenvalues()
[0.05342508803827721?, 0.779908245295057?]

returns eigenvalues in the splitting field which uses interval arithmetic. For arbitrary precision floating point numbers, there is no special implementation,

sage: matrix(RR,[[1/2,-1/4],[-1/2,1/3]]).eigenvalues()
/home/vbraun/opt/sage-5.0.beta11/local/bin/sage-ipython:1: UserWarning: Using generic algorithm for an inexact ring, which will probably give incorrect results due to numerical precision issues.
  #!/usr/bin/env python
[0.779908245295056, 0.0534250880382772]

so Sage falls back to the lapack implementation which doesn't track precision. Finally, over the symbolic ring

sage: matrix(SR,[[1/2,-1/4],[-1/2,1/3]]).eigenvalues()
[-1/12*sqrt(19) + 5/12, 1/12*sqrt(19) + 5/12]

Sage computes the symbolic answer.

Now to your question, nobody has implemented SVD for matrices with base field other than floating point numbers (the only implementation is via lapack). Probably because the real utility of the SVD is that there are efficient numerical ways to compute it, while the symbolic answer is almost always going to be a huge mess of roots.

A possible project would be to implement SVD for more base fields, for example starting with the explicit formula for 2x2 matrices.

Here's an example.

A=matrix(RDF,2,2,[1/2,-1/2,-1/2,1/2]) A.SVD()