MATLAB-problema

Messaggioda Cora_93 » 13/10/2017, 16:51

Buongiorno, sono nuovo del forum, mi trovo alle prese con un esercizio su una distribuzione di Gumbel per un'analisi sulle portate di piena.Ricavati tutti i parametri in Matlab, vorrei eseguire il plottaggio; il numero delle classi è però diverso dal numero dei campioni, il plot quindi non lo posso fare.

faj=[3 8 6 3 2 1]; (frequenza assoluta)
gj=[0.012 0.032 0.024 0.012 0.008 0.004] (densità di probabilità)
Ho pensato a questa soluzione (che funziona ma non mi piace):
g=[gj(1) gj(1) gj(1) gj(2) gj(2) gj(2) gj(2) gj(2) gj(2) gj(2) gj(2) gj(3) gj(3) gj(3) gj(3) gj(3) gj(3) gj(4) gj(4) gj(4) gj(5) gj(5) gj(6)];
praticamente faj è il numero di volte che gj deve comparire in g con il valore gj
solo che vorrei automatizzare il processo...come potri fare?
ho provato svariati cicli for ma non riesco...
so che non è una cosa complicata ma non ne esco :shock:
grazie
Cora_93
Starting Member
Starting Member
 
Messaggio: 1 di 4
Iscritto il: 13/10/2017, 16:43

Re: MATLAB-problema

Messaggioda Pazzuzu » 13/10/2017, 22:59

Di sicuro qualche informazione in più e un codice più ordinato e pulito sarebbero stati d'aiuto per chi ti volesse aiutare. Non hai nemmeno specificato chiaramente quale sia il tuo scopo. Immagino che tu debba costruire il vettore g in automatico, è corretto ? Il bello di Matlab è che si può fare tutto in pochissime righe :
Codice:
clc; clear all;

faj=[3 8 6 3 2 1];
gj=[0.012 0.032 0.024 0.012 0.008 0.004] ;

l_faj = length(faj);
g= zeros(1,sum(faj));

n=1;
for i=1:l_faj
   g(n:(n+faj(i)-1)) = gj(i);
   n=(n+faj(i));
end


Sicuramente esistono modi più eleganti per farlo, ma questo dovrebbe funzionare. Ti invito a verificare che non ci siano errori nel codice, sviste, e a commentarlo per il futuro. Attento, nella prima linea ho chiesto che venga pulito il workspace.
Pazzuzu
Average Member
Average Member
 
Messaggio: 229 di 536
Iscritto il: 08/05/2011, 16:29

Re: MATLAB-problema

Messaggioda Cora_93 » 14/10/2017, 09:57

Ringrazio per l'aiuto e mi scuso per non aver postato bene le cose, è che mi sembrava inutile postare il codice senza allegare i txt che uso, comunque il ciclo postato funziona. Io non so per quale motivo mi ero messo in testa di creare un vettore che g in maniera diversa, ma era più articolato e complesso e soprattutto non stava funzionando :roll: ...in realtà così era la prima cosa che doveva venire in mente :oops: . Grazie mille
Testo nascosto, fai click qui per vederlo
clc
clear
close all
Q=load('qq.txt');
% Q=load('Portate.txt');
m=length(Q);
Qm=sum(Q)/m;
Qsort=sort(Q);
scarti=Q-Qm;
sdv=sqrt((sum((scarti).^2))/(m-1));
M2q=(sum((scarti).^2)./(m-1));
M3q=(sum((scarti).^3)./(m-1));
M4q=(sum((scarti).^4)./(m-1));
cvq=sdv/Qm;
Gq=M3q/sdv^3;
kq=M4q/sdv^4;
nk=1+1.33*log(m); nk=ceil(nk);
Qup=88; Qdown=22; deltaQ=(Qup-Qdown)/nk;
Qclassi=Qdown:deltaQ:Qup;
Qk(1)=Qdown;

for i=2:nk
Qk(i)=Qk(i-1)+deltaQ;
end

c(1)=Qk(1)+0.5*(Qk(2)-Qk(1));
for i=2:nk
c(i)=c(i-1)+deltaQ;
end
faj=hist(Q,c);
frj=faj/m;

FRj(1)=frj(1);
for i=2:nk
FRj(i)=FRj(i-1)+frj(i);
end
gj=frj/deltaQ;

g= zeros(1,nk);
n=1;
for i=1:nk
g(n:(n+faj(i)-1)) = gj(i);
n=(n+faj(i));
end

alpha=pi/(sdv*sqrt(6));
eps=Qm-0.4501*sdv;

yi=alpha*(Qsort-eps);
for i=1:m
F(i)=i/(m+1);
end
e=exp(1);

P=e.^(-e.^-(yi));
Tr=1./(1-P);
p=(e.^(-yi)).*P*alpha;
y=-log(-log(F));
Tremp=1./(1-F);
Qi=y./alpha + eps;
TRFI=[5 10 20 50 100 200 500 1000];
PFI=1-1./TRFI;
YFI=-log(-log(PFI));
QFI=eps+(YFI./alpha);

figure(1)
hist(Q,c);
axis ([20 90 0 10])
ylabel('Frequenze assolute di classe (faj)')
xlabel('Classi')
ch= findobj(gca,'Type','patch');
set(ch,'FaceColor',[0 0.5 1],'EdgeColor','w');

figure(2)
bar(c,gj,1,'FaceColor',[0 0.8 0.4],'EdgeColor','w');
ylabel('Densità di frequenza (gj)')
xlabel('Classi')

figure(3)
plot(yi,Qsort,'m-o','Linewidth',1.5);hold on
plot(yi,Qi,'bd','Linewidth',1.5);hold off
axis ([min(yi)-1 max(yi)+1 min(Q)-5 max(Q)+5 ])

figure(4)
plot(Qsort,p,'m-','Linewidth',1.5);hold on
plot(Qsort,g,'b-','Linewidth',1.5);hold off

figure(5)
plot(Qsort,P,'m-o','Linewidth',1.5);hold on
plot(Qsort,F,'bd','Linewidth',1.5);hold off


P.s. si ho il difetto che il codice lo commento solo una volta che ho scritto la "versione definitiva più pulita".
Cora_93
Starting Member
Starting Member
 
Messaggio: 2 di 4
Iscritto il: 13/10/2017, 16:43


Torna a Ingegneria

Chi c’è in linea

Visitano il forum: Nessuno e 1 ospite