Discussion :

[jhonnyBravo]

Bonjour,

Sous matlab je voudrais afficher en 3D la fonction Bm suivante:

N = 200
mmax = 32
k(m) = N/mmax * m

Bm = sum pour p variant de k(m) à N, de nchoosek(N, p) * (m/mmax)^p * (1-m/mmax)^(N-p)

avec sur l'axe des x: m variant de 1 à 32
avec sur l'axe des y: k(m) qui prend les valeurs N/mmax, N*2/mmax, ..., N
et enfin la valeur Bm sur l'axe des z.

Voila si qqn pouvais me mettre sur la bonne voie ?

Merci

[Jerome Briot]

Pourrais-tu être plus précis ? Quel est le problème ?

As-tu déjà calculé Bm ?
Veux-tu afficher un nuage de points ou une surface ?

As-tu déjà essayé de coder quelque chose ?

[jhonnyBravo]

la fonction Bm est défini sur les entiers m=1, 2, 3, ... 32.
les valeurs k(m) se déduisent de m.
Donc si j'ai au moins l'affichage de points localisé ca me va. (apres si on peut avoir mieux du genre une nape (par interpolation donc) je suis preneur aussi.
Oui j'ai déjà codé qqch en c++ qui me permet de calculer Bm donc cest qqch qui se calcul cest ok.

Je fais appel à matlab pour afficher en 3D et en qqs lignes.

Sinon en matlab j'ai cherché comment calculer les coef binomiaux : nchoosek.

that's all

[Jerome Briot]

Dans un premier temps, utilise MESHGRID pour générer tous les couples (x,y) correspondant aux valeurs calculées de Bm et utilise plot3 pour les afficher :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
 
Bm=...;
[x,y]=meshgrid(...,...);
 
figure
plot3(x(:),y(:),Bm(:))

On verra ensuite pour améliorer l'affichage

[jhonnyBravo]

merci pour tes indications

alors : [X,Y] = meshgrid(1:1:32, 200/32:200/32:200) pour générer (x, y)
puis Y = floor(Y) car il faut des valeurs entieres.
mais comment écrire
Bm = sum pour p variant de k(m) à N, de nchoosek(N, p) * (m/mmax)^p * (1-m/mmax)^(N-p)

en matlab ?

(en fait est ce possible d'éviter l'écrire d'une boucle pour calculer les Bm, en utilisant une fonction de matlab?)

thx

[Jerome Briot]

Pour le début, tu peux faire :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
N = 200;
mmax = 32;
k=linspace(N/mmax,N,mmax);

Ensuite j'ai un doute :

Le second argument de NCHOOSEK dois être un entier. Or k (donc p) n'est pas constitué que d'entiers => bug

[Jerome Briot]

(en fait est ce possible d'éviter l'écrire d'une boucle pour calculer les Bm, en utilisant une fonction de matlab?)

Je ne pense pas (à cause de l'appel à NCHOOSEK). Mais tu peux préallouer la matrice Bm pour optimiser ton code

=> Qu'est-ce que la préallocation de mémoire ?

[jhonnyBravo]

j'ai essayé ca :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 
N=200;
mmax=32 ;
 
Bm=zeros(1, 32) ;
for m =1:32
    km = floor(N*m/mmax)  ;
    for p = km:N
        Bm(m) = Bm(m) + nchoosek(N, p) * (m/mmax)^p * (1-m/mmax)^(N-p) ;
    end
end
 
[x,y] = meshgrid(1:1:mmax, N/mmax:N/mmax:N) ;
y = floor(y)  ;
figure
plot3(x(:),y(:),Bm(:))

j'obtiens ca :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
In nchoosek at 55
  In mybin at 9
Warning: Result may not be exact. Coefficient is greater than 10^15,
         and is only good to 15 digits.

les nombres sont trop grands, en c++ j'utilisé une libraire pour faire du calcul sur les grands nombres.

Une solution ?

[Jerome Briot]

en c++ j'utilisé une libraire pour faire du calcul sur les grands nombres.

Tu devrais pouvoir englober ta librairie dans un fichier Mex : MEX-files Guide

[jhonnyBravo]

merci, je m'en suis sorti, c'est zolie

[Jerome Briot]

je m'en suis sorti

Comment ?

[jhonnyBravo]

Pour certaine valeurs matlab s'en sort très bien (jai comparé avec les résultats obtenus avec mon code c++) et à la place de nchoosek
j'ai utilisé

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 
% binomial.m by David Terr, Raytheon, 5-11-04
 
% Given nonnegative integers n and m with m<=n, compute the 
% binomial coefficient n choose m. 
 
function bc = binomial(n,m)
 
if 2*m <= n
    m1 = m;
else
    m1 = n-m;
end
 
if n <= 100
    size = abs( n-m1 ) + 1;
    ps = pascal( size );
    bc = ps( size, m1+1 );
else    % suggested by Greg von Winckel to improve speed for large input
    bc = exp( gammaln(n+1) - gammaln(n-m+1) - gammaln(m+1) );
end

téléchargeable sur le site de mathworks.

et l'algo est en fait:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
 
N=200;
mmax=32 ;
 
Km = zeros(1, 32);
for m=1:32
    Km(m) = floor(N*(m/mmax)^1); 
end
 
Bm=zeros(32, 32) ;
for m =1:32
    for i=1:32 
        for p = Km(i):N
            Bm(i, m) = Bm(i, m) + binomial(N, p) * (m/mmax)^p * (1-m/mmax)^(N-p) ;
        end
        Bm(i, m) = -log10(mmax*Bm(i, m));
    end
end
 
 
[x,y] = meshgrid(1:1:mmax, N/mmax:N/mmax:N) ;
y = floor(y)  ;
figure
plot3(x(:),y(:),Bm(:), '+')

M'enfin ce qui pourrait servir à quelquun dautre aurait été comment gérer le calcul sur les grands nombres, question ouverte.

[Jerome Briot]

Quelques modifications à ton code :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
N=200;
mmax=32 ;
 
Km=floor(linspace(N/mmax,N,mmax));
 
Bm=zeros(32, 32) ;
for m =1:32
    for i=1:32 
        for p = Km(i):N
            Bm(i, m) = Bm(i, m) + binomial(N, p) * (m/mmax)^p * (1-m/mmax)^(N-p) ;
        end
        Bm(i, m) = -log10(mmax*Bm(i, m));
    end
end
 
figure  
 
[x,y] = meshgrid(1:1:mmax, Km) ;
plot3(x(:),y(:),Bm(:), 'r+')
 
hold on
 
surf(1:1:mmax,Km,Bm)

[jhonnyBravo]

thx