This should really work but it does not :-(  frustrating!

solve(int(Mpdf, x[1] = -infinity .. xx[1], x[2] = -infinity .. xx[2]) = 0.05, {xx[1], xx[2]})

r is the Array

with(DocumentTools);
Do( %MathContainer1 = r);

For a univariate distribution we can solve for the 5-prob Quantile as follows:

restart:
with(Statistics):

int(PDF(Normal(0, 1), x), x = -infinity .. xx);
CDF(Normal(0, 1), x);

Quantile(Normal(0, 1), 0.5e-1);
solve(int(PDF(Normal(0, 1), x), x = -infinity .. xx) = 0.5e-1, xx);


But does that work for a multivariate distribution?

restart:
with(Statistics):
with(plots):
randomize():

X := proc (r) local x, t;

x[0] := .25:

for t to 200 do

x[t] := r*x[t-1]*(1-x[t-1])

end do:

[seq(x[t], t = 1 .. 200)]

end proc:

listplot([seq([r, X(r)[200]], r = 2.4 .. 4, 0.1e-2)], title = "Bifurcation Diagram", labels = ["r", "x[t]"]);

            x

I think I managed to figure it out myself.

restart:
with(combinat):
for i to 5 do x[i] := permute([`$`(1, i), `$`(-1, i)], i) end do:
x[1]; x[2]; x[3]

I think that per definition chaos is deterministic ie no-randomness-non-stochastic.
I also think that the unpredicatble nature of chaotic systems comes from the
Butterfly effect http://en.wikipedia.org/wiki/Butterfly_effect

For example for the logistic map we get:

restart:
with(Statistics):

x[0] := 0.1:

for i to 50 do
x[i] := 4*x[i-1]*(1-x[i-1])
end do:

LineChart([seq(x[i], i = 1 .. 50)]);

where the dynamics will look different depending on the value of x[0].
However, I like the way you are thinking. It is an intresting question maybe there
exist another answer?!

I have found some extra informastion here:

http://www.physicsforums.com/showthread.php?t=469212

A matrix is PSD if and only if it's principal minors are nonnegative. A 2x2 matrix
has three principal minors - the diagonal elements, and the determinant.

So x,z >= 0, and xz - y^2 >=0.

I guess xz - y^2 is the determinant which by definition has to be >=0.


Is this the purpose of the cone constraint to force the euclidean distance matrix
to be PSD??

Thanx for that Robert! In the end I did like this:

PS:=4:
seq(v[i] >= w[i], i = 1 .. NC), seq(v[i] >= -w[i], i = 1 .. NC),
add(v[i], i = 1 .. NC) <= 2*PS, add(w[i], i = 1 .. NC) = 0;

Then I get 4 long (1) and 4 short (-1) positions and the rest 0.

I think I will answer my own question:

assume = integer   seq(w[i] >= -1, i = 1 .. NC)    seq(w[i] <= 1, i = 1 .. NC)

For example:

LPSolve(add(PP[i]*w[i], i = 1 .. NC)-add(Draw[i]*w[i], i = 1 .. NC), {seq(w[i] >= -1, i = 1 .. NC), seq(w[i] <= 1, i = 1 .. NC), add(w[i], i = 1 .. NC) = PS}, maximize = true, assume = integer);

Humm quite simple! Why did I not think of that previously...hummm

@alex_01 Knowing nothing about what you are trying to do, I tried doing

eval(ob,{x[1]=0,x[2]=0,x[3]=0});

which results in  0.8571428571/x[4], which means that with the constraint con1 you can
make this as small (and negative) as you like.

Did you leave something out?

The reason is because we are allowing for long and short positions hence a short position
can completely offset a long poition which means that the problem can appear unbounded.

This is another concern; Does there exists a simple way to limit the long/short possitions ie
cardinal constraint in QP?

no I dont think so. Well, I am trying to do some portfolio optimization.
The example is taken from:  Solving Optimization Problems in Maple 9.5

http://www.maplesoft.com/applications/view.aspx?SID=4593


I am trying to expand on that example. That example minimize the portfolio variance
given a target expected return ( ob=(Transpose(X).Q.X)  con= add(X[i]*ER[i], i = 1 .. n) >= G)
I want to maximize the risk adjusted return (ob=Transpose(X).ER)/(Transpose(X).Q.X)) )
hence I need to find a way to rewrite the objective function so it remains quadratic.

n=number of stocks
C= amount of available capital
G= Target portfolio return

X = column vector with stock weights (we are trying to find optimal x[1], x[2]....etc)
ER = column vector with expected return
Q= Covariance Matrix return

con1= the sum of all our positions must be equal to or less than C
con2 = the expected portfolio return must be equal to or bigger than T

I know that you can solve the knapsack problem with dynamic programming as seen below.
But my somewhat primitive written algorithm is slow when #variable and the constraint is increased.

restart:

con := 18:
w := [4, 6, 5, 7, 3, 1, 6]:
v := [12, 10, 8, 11, 14, 7, 9]:

for j to con do
for i to nops(v) do

if `and`(i = 1, j < w[i]) then V[i, j] := 0; K[i, j] := 0
elif `and`(i = 1, j >= w[i]) then V[i, j] := v[i]; K[i, j] := 1 end if:

if `and`(i <> 1, j < w[i]) then V[i, j] := V[i-1, j]; K[i, j] := 0
elif `and`(i <> 1, j = w[i]) then V[i, j] := max(v[i], V[i-1, j])
elif `and`(i <> 1, j > w[i]) then V[i, j] := max(v[i]+V[i-1, j-w[i]], V[i-1, j]) end if:

if `and`(`and`(i <> 1, j = w[i]), v[i] > V[i-1, j]) then K[i, j] := 1
elif `and`(`and`(i <> 1, j = w[i]), v[i] < V[i-1, j]) then K[i, j] := 0 end if:

if `and`(`and`(i <> 1, j > w[i]), v[i]+V[i-1, j-w[i]] > V[i-1, j]) then K[i, j] := 1
elif `and`(`and`(i <> 1, j > w[i]), v[i]+V[i-1, j-w[i]] <= V[i-1, j]) then K[i, j] := 0 end if:

end do:
end do:

interface(rtablesize = 20):

ValueM := Matrix(nops(v), con, proc (i, j) options operator, arrow; V[i, j] end proc);
KeepM := Matrix(nops(v), con, proc (i, j) options operator, arrow; K[i, j] end proc);

cc := con:
sel := []:

for j from nops(v) by -1 to 1 do

if `and`(K[j, cc] = 1, w[j] <= cc) then sel := [op(sel), j]; cc := cc-w[j] else next end if

end do:

[[Total*value = add(v[sel[i]], i = 1 .. nops(sel))], [Total*weight = add(w[sel[i]], i = 1 .. nops(sel))], selection = sort(sel, `<`)];
 


ValueM :=

    [0  0   0  12  12  12  12  12  12  12  12  12  12  12  12  12  12  12]
    [                                                                    ]
    [0  0   0  12  12  12  12  12  12  22  22  22  22  22  22  22  22  22]
    [                                                                    ]
    [0  0   0  12  12  12  12  12  20  22  22  22  22  22  30  30  30  30]
    [                                                                    ]
    [0  0   0  12  12  12  12  12  20  22  23  23  23  23  30  31  33  33]
    [                                                                    ]
    [0  0  14  14  14  14  26  26  26  26  26  34  36  37  37  37  37  44]
    [                                                                    ]
    [7  7  14  21  21  21  26  33  33  33  33  34  41  43  44  44  44  44]
    [                                                                    ]
    [7  7  14  21  21  21  26  33  33  33  33  34  41  43  44  44  44  44]


                [0  0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1]
                [                                                    ]
                [0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1]
                [                                                    ]
                [0  0  0  0  0  0  0  0  1  0  0  0  0  0  1  1  1  1]
                [                                                    ]
 KeepM := [0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  1  1  1]
                [                                                    ]
                [0  0  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1]
                [                                                    ]
                [1  1  0  1  1  1  0  1  1  1  1  0  1  1  1  1  1  0]
                [                                                    ]
                [0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0]


     [[Total value = 44], [Total weight = 18], selection = [1, 2, 3, 5]]

I think I will answer my own question:


with(combinat):
choose([1, 2, 3, 4, 5, 6, 7]);

Finaly, I question that I think I can solve! 
It was actually alec that showed me how to do it and
now I am showing you.

restart:
f := (i, j)-> w[i]*w[j]*cov[i,j]:
C := Matrix(4, f, shape = symmetric);
Portfolio Variance = add(i, `in`(i, C));
        

    Portfolio Variance = w[1]  cov[1, 1] + 2 w[1] w[2] cov[1, 2]

                                                             2          
       + 2 w[1] w[3] cov[1, 3] + 2 w[1] w[4] cov[1, 4] + w[2]  cov[2, 2]

                                                             2          
       + 2 w[2] w[3] cov[2, 3] + 2 w[2] w[4] cov[2, 4] + w[3]  cov[3, 3]

                                     2          
       + 2 w[3] w[4] cov[3, 4] + w[4]  cov[4, 4]