We present the study of the radiative penguin Bs0 → ϕγ and Bs0 → γγ decays and the resonant B+ → K+h → K+γγ decays. We use data samples recorded at the Υ(4S) and Υ(5S) resonances with the Belle detector at KEKB, an electron-positron collider located in Tsukuba, Japan. The Υ(4S) sample corresponds to an integrated luminosity of 492 fb-1 and contains 535 million BB pairs. At the Υ(5S) resonance, we use a sample of 23.6 fb-1 containing 2.8 million Bs0 mesons. Penguin decays are loop-induced diagrams involving the heaviest particles of the Standard Model, the model describing, successfully so far, all particles and interactions, except the gravitation. Penguin decays are sensitive to physics beyond the Standard Model: particles foreseen by theories extending the Standard Model, whatever their masses are, can enter the loop and modify physics observables away from their Standard Model expectations. We report the first observation of a radiative penguin decay of the Bs0 meson in the Bs0 → ϕγ mode with a significance of 5.5 standard deviations. We measure Β(Bs0 → ϕγ) = (57 -15+18 (stat)  -11+12 (syst)) × 10-6 in agreement with the expectation of the Standard Model. We do not observe any significant Bs0 → γγ signal and we compute an upper limit at the 90% confidence level on its branching fraction of Β(Bs0 → γγ) < 8.7 × 10-6. This limit is about six times more stringent than the previously published one. However, it is still about one order of magnitude larger than the expectation of the Standard Model and still above expectations of theories beyond the Standard Model. For the resonant B+ → K+h → K+γγ decays, we search for decays where the h particle can be a η, η', ηc, ηc(2S), χc0, χc2 or a J/ψ meson, or the X(3872) particle discovered in 2003 by the Belle collaboration. We observe the modes with h = η and η'. We obtain an evidence of the mode with h = ηc; this is the first time that a B+ → K+ηc signal is seen in the K+ γγ final state. We measure Β(B+ → K+η → K+γγ) = (0.87 -0.15+0.16 (stat)  -0.07+0.10 (syst)) × 10-6, Β(B+ → K+η' → K+γγ) = (1.40 -0.15+0.16 (stat)  -0.12+0.15 (syst)) × 10-6, Β(B+ → K+ηc → K+γγ) = (0.22 -0.07+0.09 (stat)  -0.02+0.04 (syst)) × 10-6, with significances of 7.3, 13.8 and 4.1, respectively. For the other modes, we obtain limits on their branching fractions. We also measure or set limits on the branching fractions of the h → γγ decays for the modes where Β(B+ → K+h) has been measured elsewhere. We set for the first time an upper limit at the 90% confidence level on the branching fraction of the decay of the X(3872) particle into two photons of Β(X(3872) → γγ) < 1.1%.